High-energy evolution to three loops

The Balitsky-Kovchegov equation describes the high-energy growth of gauge theory scattering amplitudes as well as nonlinear saturation effects which stop it. We obtain the three-loop corrections to the equation in planar N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{N} $$\end{document} = 4 super Yang-Mills theory. Our method exploits a recently established equivalence with the physics of soft wide-angle radiation, so-called non-global logarithms, and thus yields at the same time the threeloop evolution equation for non-global logarithms. As a by-product of our analysis, we develop a Lorentz-covariant method to subtract infrared and collinear divergences in crosssection calculations in the planar limit. We compare our result in the linear regime with a recent prediction for the so-called Pomeron trajectory, and compare its collinear limit with predictions from the spectrum of twist-two operators.


Introduction
In high-energy scattering, the aspect of a particle depends on the energy scale at which it is probed. In hadronic collisions this effect can be seen in the well known energy dependence of parton distribution functions. The energy dependence can be accessed in a more detailed way by looking at less inclusive observables, for example ones probing correlations between very different rapidities, opening a window on the transverse structure of the projectile.

JHEP02(2018)058
The three-loop calculation in this paper is enabled by a recently established correspondence with the physics of wide-angle soft radiation, sometimes called "non-global logarithms." Consider the QCD decay of a color singlet state like a virtual photon or Z boson, with energy Q. A representative observable, sensitive to soft wide-angle radiation, is the probability to not find radiation with energy above a cutoff µ within an exclusion region R (see figure 1). If the cutoff µ is low, this probability is small and controlled in the planar approximation ('t Hooft limit N c → ∞) by the Banfi-Marchesini-Smye evolution equation [20]: (1.1) This resums large logarithms log Q µ . Here α ij ≡ 1−cos θ ij 2 , and the subscripts denote the angles of outgoing partons; the dipole U ij = 1 Nc Tr[U (θ i )U † (θ j )] is a function of two angles which can be interpreted (see below) as the trace of a color dipole at angles θ i and θ j .
The basic physics of this equation is that the color flow, and therefore the energy flow, is affected by radiation of an extra gluon at angle θ 0 . The observable, through the exclusion region R, is encoded by the infrared boundary condition that U ij = 0 when either i or j are in R. Qualitatively, the evolution leads to an increased effective size of the exclusion region, as radiation near the allowed boundaries become more and more in danger of leaking out. 1 This equation is mathematically equivalent to the Balitsky-Kovchegov equation, which governs the rapidity dependence of perturbative high-energy scattering near the forward direction. In this context, the trivial fixed point U ij =1 represents a transparent target, which is unstable: by linearizing in the departure (1 − U ), which gives the BFKL equation, one finds a growing solution known as the BFKL Pomeron. The nonlinear term then accounts for a class of saturation effects which stop the growth (locally in the transverse plane) toward the attractive, opaque, fixed-point U ij =0.
The nonlinear term in both equations share a similar physical origin: in both cases one is interested in the probability that something does not happen, while many possibly complicated things may happen [21]. Indeed, to describe the probability to not radiate in a certain region, one must keep track of all allowed radiation, which is what the nonlinear term of eq. (1.1) produces. Similarly, in near-forward scattering, one measures the probability for a projectile to not be destroyed at a given impact parameter. The two evolutions share other physical similarities: both are dominated by soft gluons, and both feature "opaque" and "transparent" regimes.
Given these similarities, it seems natural to expect a relationship between these two problems. The geometry is however different. To establish a rigorous map turns out to require a conformal transformation [22,23], which equates detector measurements at infinity with the physics of a fast particle crossing a Lorentz-contracted target (also known as a shockwave). This had been used notably by Hofman and Maldacena and others to describe detector measurements in conformal field theories [22,24,25] and at the same JHEP02(2018)058 time gain new insight into high-energy scattering. This conformal transformation is just the stereographic projection of a two-sphere onto the transverse impact parameter plane: Here m is an arbitrary mass scale and η is rapidity. Under this dictionary, the Banfi-Marchesini-Smye equation (1.1) becomes precisely the Balitsky-Kovchegov equation, as was noted early on [26]. In this paper we will exploit this correspondence and work exclusively on the nonglobal logarithm side, which is technically advantageous due to a body of knowledge on the infrared and collinear factorization of amplitudes and cross-sections. This correspondence was emphasized and tested explicitly at two-loops in [23], where the full two-loop BFKL/BK equation (including running coupling effects and non-planar corrections) was re-derived starting from non-global logarithm problem.
The evolution equation (1.1) at finite coupling is best viewed as a renormalization group (RG) equation: where the color density matrix, or weighted cross-section σ[U ], is defined operationally by weighing each final state parton by a color rotation U (θ i ) [23] (see also [27]). These color rotations can be understood as Wilson lines U (θ i ) accounting for the effect of more infrared radiation (these Wilson lines connect the decaying state in the matrix element and its conjugate) [28,29]. In the planar limit, the color factors reduce to products of color dipoles and the color density matrix simplifies to a single function U ij of only two angles, as shown in figure 1, which illustrates the "U U " term in eq. (1.1). In both eqs. (1.1) and (1.3), µ is an infrared cutoff below which all radiation is inclusive. In our practical calculation we will work within dimensional regularization to D = 4 − 2 ( < 0). Then the cutoff µ appears in a renormalization procedure. Following standard procedure, this is equivalent to integrating the RG equation from the deep infrared: where, writing g 2 (λ) = g 2 (µ)(λ/µ) −2 + O(g 4 ) for the running coupling in D dimensions, one can see that the integral produces 1/ poles. The subtraction then cancels the poles in the bare amplitude so as to make σ ren [U ; µ] finite as →0. That the divergences exponentiate in precisely this way was proved to all orders in ref. [23], exploiting known results on the factorization of soft partons [30,31]. The upshot of eq. (1.4) is that we can use the 1/ poles in the dimensionally regulated weighted cross-section to read off the nonglobal-logarithm/Balitsky-Kovchegov kernel K. Note that this is identical to the standard procedure to extract (ultraviolet) anomalous dimensions of local operators, by using their 1/ poles. The fact that divergences (either infrared or ultraviolet) are controlled by renormalization group equations is of course due to the Wilsonian decoupling between physics at different scales. Figure 1. Soft wide-angle radiation: radiation is allowed in some region but excluded in another.

JHEP02(2018)058
To keep track of the allowed radiation we use a color density matrix, defined by applying an angledependent color rotation U (θ i ) between the matrix element and its conjugate for each final state particle. In the planar limit this configuration reduces to a product of two color dipoles.
This paper is organized as follows. In section 2 we introduce useful notations for soft currents and phase space integrals. In section 3 we revisit the two-loop calculation, improving on previous treatments by introducing a scheme where Lorentz symmetry is manifest at each step; under the correspondence with the Regge limit, this is equivalent to maintaining the conformal symmetry of the BK equation. In section 4 we perform the threeloop calculation, paying special attention to the combinatorics of subdivergences and their cancellations, culminating in the final result for the nonlinear equation in subsection 4.6. In section 5 we analyze its linearized limit, compute its eigenvalues, and compare it with integrability predictions. Finally section 6 contains our concluding remarks. In three appendices, we record the one-loop double soft current squared (appendix A), we detail our algorithm to compute finite angular or transverse integrals (appendix B), and record the three-loop eigenvalue (appendix C).

Notations
The calculation of K requires squared matrix elements for emitting soft partons off two color-correlated parents ("dipole"), the so-called soft currents. The evolution equation, just like the soft currents, is universal and does not depend on details of the underlying short-distance process, only on the color charges and angles of the outgoing partons. Final states are then weighted, in the planar limit, by a product of color dipoles (see figure 1). For each such product, it is useful to pull out a universal factor which accounts for its dimensionality and most singular limits. We thus write the contribution from the soft current with n soft partons to the n-loop cross-section, starting from a parent dipole U 12 along directions p 1 and p 2 , as: Here loops are counted in powers of g 2 ≡ g 2 YM Nc , and phase-space integrals are normalized as p 0 ≡ 16π 2 µ 2 d 3−2 p 0 (2π) 3−2 2p 0 0 . For the Mandelstam invariants and their multi-index generalizations, we let: All invariants will always be positive (timelike), since we assume a color singlet initial state. Naturally for our setup, on-shell momenta will be split into an energy a 0 and angular parts β 0 : p µ 0 = a 0 β µ 0 where β µ 0 = (1, n 0 ) is null. The Lorentz-invariant phase space measure correspondingly splits into an energy and angular parts: For angles we write α ij = 1−cos θ ij 2 , which runs between 0 and 1. Throughout, we will use the subscripts 0, 0 , 0 to index radiated gluons.
The various factors of 2 in our definitions have been chosen to simplify limits and preclude unnecessary (log 2)'s in integrated expressions. For example, for one soft gluon, F [1 0 2] is the square of the well-known eikonal soft current. Including the factor T a T a /N c = 1/2 from the color sum, this evaluates to For two soft partons one needs the square of the double soft current, described for example in [30]. The result after squaring it and including the fermions and scalars of N = 4 SYM can be borrowed from formulas of ref. [23] (section 3), also rederived below in subsection 4.1: F [1 00 2] = 1 + s 12 s 00 + s 10 s 0 2 − s 10 s 02 2s 1(00 ) s (00 )2 . (2.5) One can easily verify that this factorizes in soft limits: Our three-loop computation builds on the one-loop corrections to F [1 00 2] and the tree-level three-parton amplitude F [1 00 0 2] , which will be efficiently obtained as described below.
3 Two-loop evolution: fixing a convenient scheme The next-to-leading order correction to the Balitsky-Kovchegov equation was obtained in QCD and N = 4 SYM in [9,32]. It was postulated [33] and verified explicitly that the same kernel governs non-global logarithms [23].
In the latter reference, soft partons were organized in terms of their energies. Because "energy" is not Lorentz invariant, this scheme did not manifest Lorentz invariance, which had to be restored manually through a finite renormalization (guaranteed to exist given the

JHEP02(2018)058
Lorentz invariance of the underlying theory). This is formally similar to the transformation used to reach to the so-called conformal scheme in the BK literature [32]. Indeed the mapping (1.2) interchanges the Lorentz and conformal symmetry of the two problems.
Here we improve on this by using explicitly Lorentz-invariant cutoffs. To fully define the scheme in which our three-loop result will apply, we thus quickly revisit the two-loop calculation.

One-loop in Lorentz-invariant form
The idea is to define the evolution so that its exponentiation generates the emission probability of a soft gluon, in the soft approximation, integrated over a complete phase space region bounded by a Lorentz-invariant cutoff. For example, at one-loop, we define the anomalous dimension K (1) so that its integral (following the first term in the expansion of eq. (1.4)) matches the emission amplitude given in eqs. (2.1a), (2.4): where θ(x < y) is a step function forcing x to be smaller than y, and Q 2 [1 0 2] ≡ s 10 s 02 s 12 defines our cutoff. From this definition one can see that Q [1 0 2] is proportional to the energy of the radiated gluon. Physically, Q [1 0 2] is the absolute value of its transverse momentum in a frame where the parents p 1 and p 2 are back to back. (This ordering variable has been used in many other contexts, see for example [34].) It is the only Lorentz invariant scale that depends on the direction but not the energies of the parent partons.
To find K (1) from the definition (3.1), we simply identify the integration over the energy component of p 0 (called a 0 in eq. (2.3)) with that over the ordering scale λ. More precisely, λ is proportional, but not equal, to the energy a 0 , because of the angle dependence of Q [1 0 2] : Inserting this change of variable into the right-hand-side of eq. (3.1) using the measure (2.3), and stripping off dλ λ λ −2 on both sides, we thus get: This of course reproduces the one-loop Banfi-Marchesini-Smye equation recorded in (1.1), except for the in the exponent, which arose because of the angular dependence of the ordering variable Q [1 0 2] . This exponent ensures exact Lorentz invariance in any dimension, not only in the → 0 limit 2 , which is critical to ensure Lorentz invariance of the higher-loop corrections to K [23]. We briefly comment on the inclusion of virtual corrections, which simply add the (−2U 12 ) term to eq. (3.1). This form is determined by the Kinoshita-Lee-Nauenberg (KLN) JHEP02(2018)058 theorem [35,36], which states that there can be no infrared divergences in a fully inclusive cross-section. This implies, in particular, that U ij = 1 is a fixed point of the evolution. At any loop order this can (and will) be used to obtain the coefficient of U 12 from that of other color structures.

Lorentz invariant slicing of multi-particle phase spaces
To move on to higher loops, we define, similarly, scales for multiple emissions: (3.4) Similar combinations appeared already in the integration measures in eqs. (2.1). The exponents may appear unwieldy, but in practice these definitions will be very convenient because the scales of complicated processes are equal to appropriate geometric means of subprocess scales, for example: As an organizing principle, when writing the higher-loop contributions to the evolution kernel K, we make sure to completely cover the multi-parton phase space up to a cutoff in Q. Let us consider for illustration a term arising from an -loop virtual correction to the emission of two real partons ( = 0 being the relevant case for the two-loop kernel to be detailed shortly). If I denotes the corresponding soft current, the following expression integrates it over all the phase space with Q [1 00 2] below the cutoff: Importantly, the integrand will always be homogeneous, due physically to the fact that the Yang-Mills coupling is dimensionless. More precisely, within dimensional regularization, the -loop correction to the two-parton emission has an overall dimension determined by the running coupling g 2 (λ) ≈ λ −2 , raised to the power ( + 2). We thus change variable from the two energies (a 0 , a 0 ) to the overall scale λ ≡ Q [1 00 2] and relative energy τ = a 0 /a 0 . Dimensional reasoning then implies that, after factoring out the running coupling evaluated at that scale Q [1 00 2] , the integrand becomes homogeneous and depends only on the ratio τ , but not λ: With this change of variable the two-particle phase-space then factors as The integral over the scale λ precisely matches what appears in the integrated renormalization group equation (see eq. (1.4)), so by simply stripping it off we get the contribution to the kernel K ( +2) , simply generalizing eq. (3.3). Note, importantly, that the equality is JHEP02(2018)058 exact to all orders in and holds not only for the 1/ poles. The only assumption is that the integrand I is computed in the leading soft approximation where the energy scales of the parent partons does not enter, e.g. I is the standard soft current comes from soft currents, The generalization to more partons is immediate: in our slicing scheme we will always get integrations that depend only over relative energies τ , with an -free measure.
(Strictly speaking, the identity (3.8) is only valid when the integrand I vanishes in its factorization limits τ → 0, ∞ so the τ -integral converges, which holds for the subtracted integrand F sub to be defined shortly.) This equivalence between scale integrals and energy integrals can also be applied in the reverse direction, to subtract the iteration of the lower-loop kernels generated by the path-ordered exponential (1.4). For example the product of two K (1) 's corresponding to the successive emission of parton 0 between 1 and 2, followed by parton 0 between 0 and 2, can be written as (3.9) where r [1 0 2] = (α 12 /(α 10 α 02 )) 1− is the angular measure in eq. (3.3).

Quick rederivation of two-loop evolution
With this technology it is now rather straightforward to re-derive the two-loop evolution equation. Let us start with the contribution from two real partons. Matching with eq. (1.4), this requires the squared matrix element for two partons, minus the iteration of one-loop subprocesses. This later subtraction will neatly remove all subdivergences. There are two possible one-loop subprocesses: either p 0 or p 0 can be radiated first. The relation (3.9) allows to subtract these directly at the integrand level, by defining a subtracted soft current: Multiplying with the product of dipoles U 10 U 00 U 0 2 in eq. (2.1b) and removing the overall scale integral using eq. (3.8), we get K as a convergent integral: where the omitted terms involve virtual corrections (involving products of fewer than three U dipoles). The τ integral converges absolutely in both the τ → 0 and τ → ∞ limits thanks to the factorization of the soft current F noted in eq. (2.6). Note that we have omitted the µ upper cutoff in the step functions in F sub [1 00 2] . This is because all the 1/ poles at two loops come from the infrared region where p 0 ∼ p 0 µ, where this cutoff plays no role [23]. The region near the upper cutoff only affects the twoloop amplitude by a finite amount, thus affecting the evolution starting only from three loops (in a way which can be systematically accounted for, see eq. (4.21c)). The τ -integral in (3.11), using the explicit expression (2.5), involves only elementary integrals and gives JHEP02(2018)058 a simple angular function α 12 α 00 α 10 α 02 + 1 + α 12 α 00 α 10 α 0 2 − α 10 α 02 log α 10 α 0 2 α 10 α 02 .
(3.12) Finally we turn to the virtual corrections, which can have color factors U 10 U 02 or U 12 . They are strongly constrained by physical principles: Lorentz invariance, the absence of collinear singularities, and the KLN theorem. A simple way to solve these constraints is to add (U 10 U 02 + U 10 U 0 2 ) to the color factor in eq. (3.11), which automatically removes collinear singularities when 0 0 (where U 00 → 1) and fulfills KLN. By Lorentz invariance, the remainder is then determined up to a single multiple of one-loop: (3.13) The coefficient γ (2) K can be fixed by matching a certain limit controlled by the cusp anomalous dimension (see section 5): γ The full two-loop planar evolution is then given as (3.13) which agrees completely with the existing result for the Balitsky-Kovchegov equation [32].

More on virtual corrections
Although they were not strictly needed to obtain the two-loop result (3.13) (having taken the two-loop cusp anomalous dimension as a known input), it is instructive to explicitly compute the virtual corrections. Learning to handle them expediently will prevent them from becoming the bane of our existence at higher loops.
Morally, the coefficient γ K is related to the one-loop correction to the single soft current, which has been obtained long ago (see [37,38]): (3.14) Here c Γ = Γ(1+ )Γ(1− ) 2 Γ(1−2 )(4π) − is a ubiquitous loop factor. This formula does not depend on the matter content of the theory. The "bare" superscript indicates that we have performed ultraviolet renormalization but have not yet subtracted the infrared divergence, to which we now turn.
Obviously this result is divergent, whereas we're trying to compute the finite coefficient γ Of course, what happens as usual is that the physics cannot depend on such a "bare" quantity but only on renormalized ones. A useful intuition here is that infrared divergences in bare amplitudes reflect that scattering states are defined in the deep infrared, and one must always use the renormalization group to evolve the amplitude back to the physical scale µ of interest, as detailed in eq. (4.9) below. This will remove all remaining 1/ poles. In the present case, the precise renormalization to use, including finite factors, follows from the other virtual contributions already included in K. First there is the U 12 term in K (1) ,

JHEP02(2018)058
predicted above using the KLN theorem, which can multiply the real part of K (1) iterated using relation (3.9): Second there are the U U terms in the two-loop ansatz (3.13): It is important to note that both these contributions are expressed in terms of the phase space of two real partons 0 and v, whereas the one-loop virtual correction to the soft current (3.14) is to be integrated over the phase space of a single partonp 0 . We thus have to match these phase spaces somehow. The crucial requirement is that the collinear singularities match at the integrand level. This requires that, in the limit where 0 and v are collinear, their total energy matches that in the virtual calculation: a 0 + a v =ã 0 . The simplest way to do this, while respecting Lorentz invariance away from the collinear limit, is to keep the angles the same,p 0 ∝ p 0 , but use Q [1 0 2] and Q [1 v 2] to define Lorentzcovariant energies for the two daughters 0 and v. Thus we match the above two corrections with the virtual one at total momentump 0 ≡ p 0 . Let us denote as f split (p 0 , p v ) the sum over the five terms in (3.15)-(3.16), or more generally any homogeneous function of p 0 , p v . After changing variable from p 0 top 0 the two-parton phase space factorizes as: , H parent (p 0 ) is an arbitrary test function, and x and 1 − x represent the (covariant) energy fractions of the two daughters.
The splitting function f split defined by the sum of (3.15)-(3.16) contains complicated angle-dependent step functions, which come both from the former equation and from those in F sub , explicited in eq. (3.10). Conveniently, up to a part that is antisymmetric in x → 1 − x and therefore cancel upon integration, all the step functions cancel out except those proportional to θ(Q [1 v 2] < Q [1 0 2] ). Keeping only these surviving terms, and decomposing the sum into two pieces for later convenience, we thus write Stripping off the integral over the radiated gluon momentump 0 in eq. (3.17), we then get the total effective soft current:

JHEP02(2018)058
with δ the integral over the splitting function: Note that, although it is defined as a complicated looking integral, δ (1) is just a constant: this had to be the case since the integral is manifestly Lorentz-invariant and an homogeneous function of three null vectors, and all such invariants are constant. Adding it to the bare matrix element (3.14) according to (3.19) then gives: in perfect agreement with the two-loop cusp anomalous dimension recorded below eq. (3.13).

Three-loop evolution
We now proceed to derive and assemble the ingredients for three-loop infrared divergences. The chief conceptual issue is to organize the subtraction of subdivergences, of which there are plenty at three loops. We would like to (and will) obtain the evolution kernel K (3) as a sum of absolutely convergent integrals involving physical building blocks (the so-called remainder function) in which we can set = 0 directly.

First ingredient: triple-soft current
The first building block is the square of the tree-level soft current for emission of three partons. This needs to be summed over all produced parton species: gluons, fermions or scalars.
The easiest way to obtain it is from the soft limit of the planar four particle integrand, which is amply documented in the literature. We fix two external legs to be in the matrix element and two in the conjugate, and sum over all (Cutkoski) unitarity cuts which separate them. In the relevant limit, where the cut internal propagators become soft, the integrand from the outer square factors out and the outermost cut propagators act as the parent dipole U 12 .
As an illustration, consider the two-loop integrand, which in planar N = 4 is a sum of two double-boxes. With the momenta labelled as in figure 2a: . (4.1) Taking p 1 , p 2 and p 0 to be on-shell with p 0 soft, this simplifies to where the first factor is recognized as just the cut of the one-loop amplitude (a scalar box).
Dividing it out leaves the dipole radiator (2.4), as expected. The other three-particle cut of the same diagram (where the cut runs south-east) adds the correct factor 2, and the rotated double-box is subleading in the soft limit.
This is in perfect agreement with the direct calculation recorded in eq. (2.5).
Having thus validated the method, it is a simple exercise to extract the square of the triple-soft current from the known 4-loop integrand. We found the 7-loop package [40] (recently extended to 8 loops [41]) particularly useful for this. To most usefully record the result, we note that its soft limits are easily predicted. There are five independent soft limits, where (by factorization) it must reduce to double-soft currents: As a cross-check, we have reproduced numerically the squared soft current (4.4)-(4.5) by a direct Feynman diagram calculation, summing up the gluon, fermion and scalar contributions, and also using the computer package [42]. For convenience, this formula, and others in this paper, is included in computer-readable format in the ancillary text file formulas.txt, attached to the arXiv submission of this paper.

Second ingredient: double-soft current and the remainder function
To obtain the one-loop correction to the double soft current in the simplest way, we take the limit of two soft partons in the known one-loop six-point amplitude. These soft partons can be of any species (gluons, fermions and scalars). Consider for example the case when the two soft gluons have the same helicity. In this case we use the one-loop correction to the MHV amplitude (four positive and two negative helicity gluons), divided by the tree amplitude [43]:  where the operation C is a cyclic rotation by one. The one-loop soft current is obtained by taking the limit where partons 2 and 3 become the soft partons 0 and 0 , and subtracting the one-loop correction to the parent four-point amplitude. In this limit, the two coloradjacent partons 1 and 4 define the parent dipole, and the other two decouple, thus giving us the soft current It is important to note that since all invariants are positive (timelike), the Feynman prescription adds an imaginary part to all logarithms: log(−s ij ) = log |s ij | − iπ.
For soft gluons of opposite helicity, as well as for soft fermions and scalars, one needs the NMHV (super)amplitude [44,45]. It may be amusing to note that the two fermions soft current is the same in QCD and N = 4 SYM, since the contributing diagrams are the same. Thus some effective supersymmetry can also be used at one loop in QCD as well.
The component formulas are somewhat involved, and in the N = 4 theory further simplifications occur when summing over particle species in the interference with the tree amplitude. For this reason, here we record only the final result of the helicity sum, e.g. the one-loop correction to the squared soft current, in appendix in eq. (A.1): We used the package in [46] to cross-check our expressions. Importantly, as was the case at two loops (and for the MHV example above), this one-loop correction is infrared divergent, while we expect the physics to depend only on renormalized, finite quantities. The standard, MS way to renormalize is to remove the JHEP02(2018)058 integral of the infrared anomalous dimension: where at one-loop γ IR = g 2 YM Nc 8π 2 log |s 10 s 00 s 0 2 | |s 12 |µ 4 for the soft current squared. This is the conventional definition of so-called hard matrix elements in the SCET literature. Although a good starting point, this is however not very convenient for us, because we would like to subtract something which has a simple representation as a phase space integral.
The governing physical principle is that the subtraction should match all singularities of the triple-real emission at the integrand level, in all (single-) soft and collinear limits. This ensures that when we add it back later all divergences will cancel cleanly pre-integration. Furthermore we would like a simple analytic form for the integrated subtraction. This can be achieved by defining Lorentz-invariant functions of three angles, like we did in section 3.4, since these automatically integrate to constants.
Let us thus consider the general problem of renormalizing an amplitude F [1 23... n] with (n − 2) soft partons. We want to renormalize it by adding, say at one-loop, a phase space integral with one additional real parton v: There are two constraints. In the limit where v is soft, the integrand should reduce to minus the square of the soft current, s iv s vj , emitted from all possible regions, minus that from the parent dipole. In collinear limits there is a similar factorization, but with the important distinction that the parent amplitude must be evaluated with the total momentum (here j = 2, . . . , n − 1 and i = j − 1, k = j + 1): (4.11) Note that the labels i and k decouple in the collinear limit. The fact that the argument of the amplitude is shifted to (p j +p v ) is the main complication since it precludes a simple multiplicative solution. To solve it, recycling the ingredients in the two-loop subtraction (3.18), we define a three-index operator Γ [i j k] which rescales p j in whatever it multiplies: Then, it is easy to see that all constraints are simultaneously solved by: Indeed, the three-index Γ v [i j k] only has collinear singularities in one region p v p j , where the spectator labels i and k decouple, so the collinear limits work out. In the soft limit it JHEP02(2018)058 and using telescopic cancellations one can also see that the first of (4.11) is fulfilled.
Using the change of variable (3.17) and the integral (3.20), the renormalization defined at one-loop by eq. (4.10) can be rewritten in a more suggestive form: where represents the real-emission correction to the soft current squared between legs i and k. We recall that δ (1) ≈ 2 2 (3.20) starts with a double pole, and g 2 (λ) ≈ (λ/µ) −2 g 2 (µ) is the D-dimensional running coupling. Physically, the renormalized amplitude is thus obtained by including amplitudes for a sequence of splittings, each with the coupling evaluated at its natural scale. (Either of the sequences in the square brackets would work, but we chose to include both and multiply the exponent by 1 2 for symmetry reasons.) The renormalized amplitude F ren is finite for any number of points. At one-loop F (1)ren is obtained from the bare result (A.1) by the simple substitution given in eq. (A.2). It turns out that F ren is closely related to another canonical finite function in N = 4 SYM: the Bern-Dixon-Smirnov remainder function [47]. This is defined by dividing the amplitude by an ansatz A BDS n+2 (essentially an exponential of the one-loop MHV amplitude), which makes it finite and dual-conformal invariant and trivializes its collinear limits. The ansatz has four parameters: three are essentially the constant, order and 2 terms in the function called f ( ) in [47], which multiplies the one-loop amplitude in the exponent, while the fourth adds a common multiplicative factor to all n-point amplitudes and cancels out for the soft current. Thus three parameters affect the soft current; comparing with eq. (4.15) it is easy to see that these three parameters are in one-to-one correspondence with the double-pole, single-pole and constant term in δ (1) , and that the infrared divergent parts match. The 0 term in our δ is slightly different because for five-partons our scheme automatically yields the cusp anomalous dimension F [1 0 2] = γ K (see section 5 for the higher-loop explanation) whereas by definition the BDS remainder is unity for five-points. Thus, with a somewhat schematic notation, we can express our renormalized soft current directly in terms of the BDS remainder to all loop orders: where R n+2 = A n+2 /A BDS n+2 is the BDS remainder (e.g. the amplitude A n+2 divided by the BDS ansatz A BDS n+2 ) and f n an explicitly known (and finite) function equal to the real part of the difference between the exponent in eq. (4.14) and the squared one-loop MHV soft current (given for n = 4 in (4.7)), which was exponentiated by BDS. As we will prove shortly, only the → 0 limit of the finite F ren , and thus also BDS remainder, is needed to get the evolution kernel K.
Finally, it is interesting to look at eq. (4.16) in the other direction, going from the soft current to the remainder. We can count the number of variables on which the soft current F ren depends for m soft partons. Each on-shell parton gives 3m degrees of freedom while JHEP02(2018)058 invariance under two Lorentz generators and one independent rescaling of the β i remove 3, giving 3(m − 1) invariants. Due to dual conformal symmetry, the k-point BDS remainder A k /A BDS k in eq. (4.16) depends on 3(k − 5) invariants, which is equal since k=m + 4. These two numbers agree! That is, dual conformal symmetry implies that the soft limit is not lossy, and we conclude that, through eq. (4.16), the → 0 limit of BDS remainder in planar N = 4 SYM uniquely determines the renormalized soft current, and vice-versa.

Nested subtractions for virtual contribution
Given the renormalized amplitude, it is natural to integrate it over relative energies to obtain a contribution to K, with suitable subtractions as was done in section 3.3: For the subtracted integrand F ren,sub it would be tempting to use again eq. (3.10), but one needs to be more careful and pay due attention to the renormalization scales of the various objects. Indeed, as is clear from the renormalization group equation (1.4), all couplings in the subtractions get evaluated at their private scales Q [i j k] , which are distinct from the common overall scale Q [1 00 2] that we assign to K ren [1 00 2] . In addition, the finite parts of the renormalization (4.14) do not match. The correct loop-level definition, which accounts for all these effects, is rather For the other subprocess we get the same but with −κ (1) . Specializing to what we need at three loops, extracting the coefficient of g 2 (Q [1 00 2] ) in F ren,sub and using that F

JHEP02(2018)058
The critical conceptual point here is that we won't need the O( ) terms in this expression. This is because the combination in eq. (4.18), in which all objects are defined to all orders in , is precisely the one which vanishes to all order in near the endpoints τ → 0 and τ → ∞ (this follows from the factorization properties of the bare amplitudes F bare ). This precludes any / effect. The extension to higher loops is clear: one just includes more terms in the expansion of δ. Also we expect only minor changes in the presence of a nontrivial β-function as in full QCD, where g 2 (λ) will now be a series in g 2 (Q [1 00 2] ).

Nested subtractions for triple real contribution
We now turn to the fully real contribution to K (3) , which is given by the IR divergent part of triple-real emission, minus the subdivergences associated with iterations of K (1) and K (2) . The basic idea is to write the subtractions as phase space integrals with step functions, exploiting (3.9) and its higher-multiplicity generalizations. In this way all energy sub-divergences (with fixed angles, as appropriate since the angles are fixed by the color rotations U ) will cancel under the integration sign. To write the result concisely, we recursively define subtracted integrands F sub , generalizing eq. (3.10). Introducing the abbreviations these are defined as: The structure is straightforward: there is one subtraction for each possible subprocess (consistent with the planar structure), and the unsubtracted F 's are given in eq. (2.5) and (4.4). Intuitively, the F sub 's are a device to compute the logarithm of F : the preceding equations can be generated (and generalized to all orders) by formally solving the equation Pe F sub = F , order by order in the number of emitted partons. As shown in section 3, what is relevant for the evolution is the integral over relative energies: Thanks to the pattern of subtractions, and to the factorization of soft currents (see eqs. (2.6) and (4.3)), F sub [1 00 0 2] vanishes in all soft limits and its energy integral at fixed angles is absolutely convergent at all orders in . One might worry that the step functions make it tricky to integrate in practice, but in fact they always multiply trivial measures like dτ /τ . Furthermore, the explicit expression (4.4) naturally splits into several individually JHEP02(2018)058 convergent pieces. For example, the piece F safe doesn't contain any step function and converges by itself. The pieces from the "1" in F [1 0 2] , F [1 00 2] and F [1 00 0 2] contain multiple step functions, but all share the trivial measure dτ /τ dτ /τ and so immediately integrate to logarithms. Finally, the five nontrivial subtractions in (4.4) naturally combine with the remaining terms in (4.21c), to produce five individually convergent integrals.
So our problem is reduced to computing finite energy integrals; these produce functions of transcendental weight 2. A good, systematic way to compute such integrals is the differential equation method described B. 3 The most difficult integrals are contained within F safe . One of them, in particular, coming from the first line below eq. (4.5), cannot be written simply in terms of the angular distances α ij , but requires associated spinors (β αβ i ≡ λ α iλβ ): Here we have used a commonly used notation for the Lorentz-invariant spinor products: i j = αβ λ α i λ β j and [i j] = αβλα iλβ j with antisymmetric. (Under the stereographic projection (1.2), these map respectively to: i j = (z i − z j ) and [i j] = (z i −z j ).) The other integrals are more elementary and produce at most dilogarithms of cross-ratios of α's.
Then the triple-real integral gives where f 1 is the special function in eq. (4.23), P exchanges labels (1, 0) and (2, 0 ) and acts

Nested subtractions for renormalization counter-terms
The final ingredient is the "add" part from the "add and subtract" game that led to the renormalized amplitude (4.18). These also have three angular integrations but one fewer color dipole U . In addition there are similar pieces inherited from lower loops, for example from the term with just U 10 U 0 2 in the two-loop evolution. It is useful to devise a notation for such terms, like G {1 00 2} , wherein the underlined index represents the angle from which a Wilson line is omitted and the curly bracket highlights the presence of a virtual parton. This is why we've split the virtual correction into two terms (G {1 v0 2} +G {1 0v 2} ) in eq. (3.18), because these two end up with different color structures and so get exponentiated at different scales: Similarly, to exponentiate the three-loop kernel (as would be needed for a putative four-loop calculation), we would need to specify where v fits within the color structures of Γ v [1 00 2] , which determines the relevant scale Q. Thus although not strictly necessary here, it is useful to account for that information because it helps show the internal logic. Thus we organize the "add" terms into three color structures: Here we only show the terms in K (3) with three angular integrations, two has been dealt with in subsection 4.3 and one will be dealt with shortly. The underlined index shows the variable whose Wilson line and radiator factor are omitted. The angular functions are the integrals over relative energies of corresponding G sub 's, The G sub 's contain two ingredients. First, there is the difference between the renormalized F ren,sub in eq. (4.18) and the corresponding bare expression: (4.28)

JHEP02(2018)058
Second, there is the subtraction of everything inherited from virtual corrections at lower loops: from the U 12 term at one-loop (3.3), and from the (U 10 U 02 + U 10 U 0 2 ) part of twoloop (3.13). To allow their subsequent exponentiation, the result is to decomposed into 3 color structures: A simple systematic color decomposition for the subtractions can be done as follows. Whenever, in a subprocess, both indices adjacent to v are the same as in the considered G sub , we weight this contribution by 1; when only one index is shared, we weight by 1 2 , and when none is shared, we weight by 0. For example, consider the following term coming from the real part of K (1) times the virtual part of K (2) : We place half of this term into G sub {1 v00 2} and half into G sub {1 0v0 2} , because v occurs between 1 and 0 . Using these rules to generate the subtractions recursively, using the same notation as in eq. (4.21c) (writing {a b . . . c} ≡ G sub {a b... d} and inserting a step function between each bracket, either curly or square), then gives where P is the parity (10)↔(0 2). This looks messy, but the upshot is that the internal logic is straightforward and the terms can be automatically generated to any desired order. (Formally, the terms can be generated by series-expanding the schematic formula

JHEP02(2018)058
Pe (F sub +G sub ) = (F + G)e K| U 12 . 4 ) This generalizes the subtractions used at two loops: the energy integral of G sub {1 v0 2} matches the U 10 U 0 2 term in eq. (3.13). Although we haven't defined the individual Γ {1 v00 2} (only their sum Γ v [1 00 2] ) we expect that a definition exists which will make the energy integrals converge absolutely for each of the color structure G sub {1 v00 2} , as this is certainly the case for the sum which is all we need here at three loops. Furthermore, by construction, the collinear singularities of K (3)add cancel exactly those of K (3) , to all orders in , so it is apparent that O( ) corrections to any kernel are not needed.
Thus we only need to compute the finite integral (4.26) with = 0. The integrated result turns out to be somewhat inelegant, so we decided to replace it by a simpler counterterm with the same collinear singularity. From inspection of the triple-real result (4.24), we find divergences as 0 0 or 0 0 , and also in the double scaling limit 0 0 0 , but not when one or two partons become collinear to 1 or 2. A simple counter-term which removes the divergence as 0 0 is: (4.32) To construct an integral that is also absolutely convergent in double collinear limits, we can easily play with the color structures, exploiting that U ij → 1 when i j. Arranging for each color factor to separately fulfill the KLN theorem (vanishing when U ij = 1), the full three-loop evolution is then written as (4.34c) below, where the difference compared to subsection 4.3 is simply: with P the symmetry (10)↔ (20 ). This is again an absolutely convergent integral which can be done at = 0, using the methods of appendix B. We find a surprisingly compact result: − 1 + α 12 α 00 α 10 α 0 2 − α 10 α 02 log 2 α 12 α 00 α 10 α 02 + 4ζ 2 log α 10 α 0 2 α 10 α 02 − 11 6 log 3 α 12 α 00 α 10 α 02 .  31)), in such a way that all cancellations are manifest at the integrand level and valid to all orders in . This allowed us to set = 0 directly in all integrals and be completely certain that we did not miss any / effect.
We have used the squared amplitude for triple-real emission and also the one-loop correction to double-real emission (related to the one-loop six-point remainder function). In addition K receives contribution from single-real emission at two-loops, and fully virtual corrections. However, it is not necessary to explicitly compute them. As mentioned already, fully virtual corrections follow simply from the KLN theorem. And by Lorentz symmetry (kept manifest at all stages of our calculation) the single-real emissions can only produce a constant γ (3) K time one-loop. As argued (and tested) in the next section, provided that the U 12 color structure appears nowhere else in our expression, what multiplies one-loop must be the cusp anomalous dimension (known to all loops [48]): Thus our final result for the three-loop BK equation, recalling the lower loop results, is: where P is the parity (10)↔(20 ), α ij ≡ |z i −z j | 2 are transverse distances and β 0 ≡ d 2 z 0 π . (Equivalently, for the non-global-logarithmic problem, the stereographic projection (1.2) gives α ij ≡ 1−cos θ ij 2 and β 0 ≡ d 2 Ω 0 4π ). The two-loop transverse function K (2) [1 00 2] was given in eq. (3.12), and the triple-real function K the effective single-virtual kernel (the sum of eqs. (4.17) and (4.33)) is given as

JHEP02(2018)058
For convenience, these formulas are reproduced in computer-readable format in the ancillary text file formulas.txt, attached to the arXiv submission of this paper. We note that eq. (4.36) is a single-valued combination of polylogarithms. That is, it does not have any branch cut for physical angles (where x andx are complex conjugate of each other:x = x * , as is easily verified). This has to be the case since the kernel represents a physical probability for radiation and there can't be multiple answers for a given set of angles. Concretely, although this is not manifest, one can verify that the series expansion of the last line around x =x = 1 contains only single-valued logarithms of the type log(1 − x)(1 −x), but log(1 − x) never appears separately from log(1 −x).

Linearized evolution and BFKL Pomeron trajectory
In many applications to the high-energy limit, especially those involving dilute targets and projectiles, the Wilson lines remain close to unity and the physics is governed by a linearized version of eqs. (4.34). Then we set and treat U ij as a small quantity. Generically, in the 't Hooft large N c limit, U ij ∼ 1/N 2 c when scattering objects made of a fixed number of partons, or for example a four-point correlator of single-trace operators. The resulting linear equation is referred to as the BFKL equation and its eigenvalue j = 1 − K is the Pomeron Regge trajectory. With this application in mind, in this section we will use the language of transverse plane and conformal symmetry, instead of the stereographically equivalent language of angles and Lorentz symmetry. Linearizing the color structures in the three loop result (4.34c) produces many terms, but these turn out to organize simply into the combination which appears already at two loops: This is due to an exact symmetry: the large N c theory is invariant under the local gauge transformations U ij → U ij e α i −α j , representing independent U(1) gauge transformations in the past and future. (Beyond the planar limit, only the global SU(N c ) past ×SU(N c ) future survives as a symmetry of the Balitsky-JIMWLK equation.) The combination (5.2) is the only one invariant under the linear transformation U ij → U ij + α i − α j , which does not contain U 12 and is invariant under the parity (10)↔ (20 ). That parity is automatic for any conformally-invariant function of four transverse points 1, 0, 0 , 2, and so not really an assumption.
The first four terms on the right of eq. (5.2) naively integrate to zero,

JHEP02(2018)058
because by conformal symmetry the z 0 integral can only produce a constant and thus cannot be antisymmetric in 1 and 2. In the first equality we have used the just-mentioned parity symmetry to trade (0↔0 ) for (1↔2). A subtlety however is that in this rewriting the middle integral fails to be absolutely convergent in the double scaling limit z 0 , z 0 → z 2 , even though the left-hand side is. Due to this, the conformal symmetry argument breaks down for z 0 = z 2 , enabling a contact term δ 2 (z 0 −z 2 ) to appear. (See appendix E of ref. [32] for explicit examples.) Taking into account the possibility of such contact terms by adding a constant C (L) , the linearization at L-loops takes the general form where at two loops K (2)lin [1 00 2] and at three loops [1 00 2] + 2 This integral, like others in this paper, is absolutely convergent. The factor of two accounts for the contribution with 0 and 0 interchanged, which produces the same result due to the parity symmetry. We computed this integral using the differential equation method explained in appendix B. The resulting function of the cross-ratios x,x (defined in eq. (4.35)) has 5 letters in its symbol: . At transcendental weight 3 there exists rather few such functions that are real and single-valued in the physical regionx = x * , in the sense explained below (4.36). We have found only three nontrivial ones O 1,2,3 . Since there is limited information content in these functions themselves, we record them in appendix in eq. (B.10) and here record the concise coordinate space expression for the BFKL kernel (u = xx, v = (1 − x)(1 −x)): Now the linear equation (5.4) can be diagonalized explicitly because its eigenfunctions are determined by conformal symmetry [49]. The eigenvalue depends on two quantum numbers: a scaling dimension ν and an (integer) angular momentum m. It can be extracted by looking at the translation invariant wavefunctions 5  A special eigenfunction is U (z i − z j ) = z i − z j , corresponding to m = 1 and ν = 0, which is a generator of the aforementioned U(1) gauge symmetry at large N c . The eigenvalue of K = 1−j must thus vanish in this case, which leads to an exact prediction for the intercept of the Odderon trajectory at large N c [19,50]: Since this property was manifest in the original starting point, i.e. eq. (5.2) before eq. (5.3) was used, in practice we will use this property to fix the constant C (L) . Translationinvariance of the trial wavefunctions enables one transverse integral to be done explicitly. This simplifies the evolution (5.4) to: where, labelling four points as {z 1 , z 0 , z 0 , z 2 } = {1, z, z − y, 0}, the translation-invariant kernel is Plugging in the wavefunction (5.7), we see that the Pomeron trajectory is the Mellin transform of H (L) (y). The parity symmetry of K (L)lin makes |y|H (L) (y) invariant under the inversion y → 1/y. The eigenvalue can thus be written as the sum of two terms, analytic in the lower-and upper-half ν-planes respectively, representing the contributions from |y| < 1 and |y| > 1. Following a common notation in the literature, we thus write the Regge trajectory j = 1 − K as:

JHEP02(2018)058
where, for L > 1, In summary, the Pomeron trajectory j(m, ν), defined conceptually as the eigenvalue of the evolution (5.4) on the eigenfunctions (5.7), is obtained at three loops by computing the Mellin transform (5.12) of the translation invariant projection (5.10) of the coordinate space kernel (5.6).

Result for the eigenvalue
To efficiently integrate (5.10) we used the differential equation method, wherein derivatives are iteratively computed and simplified using integration-by-parts identities. This method has a long and successful history in the context of dimensionally regulated Feynman integrals [51,52]. We used a variant that exploits simplifications occuring for absolutely convergent integrals, based on ideas in refs. [53,54]. Our procedure is illustrated in appendix B in a few examples. The result is an expression for H (3) (y) in terms of iterated integrals starting from the origin y = 0. In principle these iterated integrals could be rewritten in terms of polylogarithms, but we found this neither illuminating nor useful in practice. Rather, to extract the eigenvalue, we found it more efficient to perform the angular integration directly at the level of the iterated integral, using again the differential equation method to obtain the result as a iterated integral in the radial variable x = |y| 2 . The radial functions then turned out to be conventional harmonic polylogarithms. At one-and two-loop this procedure gives expressions that are very uniform for all transverse angular momentum m (m ≥ 0): Here the 'reg.' notation is an instruction to subtract all the negative powers of x (and powers of log(x) they multiply) from the series expansion of the bracket around x = 0. The harmonic polylogarithms (with omitted argument x) are defined recursively as [55,56] 14) The concise expression (5.13) for the two-loop eigenvalue is apparently new, but we have verified that it agrees, for all values of |m|, with the known result in N = 4 SYM [10,11,57]. Equations (5.13) takes the form of a Mellin transform over harmonic polylogarithms, which is well-known to give harmonic sums (see appendix C), which in the case of eq. (5.13) would have argument 1+|m|+iν 2 and 1−|m|+iν 2 . However, it is important to note the "reg"

JHEP02(2018)058
subscript in that equation, which implies that a number of powers of x terms, which grows with |m|, have to be subtracted. It would be interesting to see if the result can be usefully written as some kind of "regulated" harmonic sum.
At three-loops, although the "reg." notation seems to help, we did not succeed in finding a compact formula accounting for the full m dependence, and so here we restrict our attention to individual values of m, for example The Mellin integral in eq. (5.15) gives a practical and efficient way to compute the eigenvalue numerically for any desired value of ν. The result of the integral can also be formally expressed in terms of harmonic sums (see eq. (C.3)), although evaluating these sums for complex ν then requires an analytic continuation. In appendix C we also provide harmonic sums expressions for m = 1.
A Mathematica notebook trajectories 3loop.nb attached to the arXiv submission article allows to evaluate the eigenvalues for any m and ν.
. This is related to the square-root containing entries of the symbol of H(y) recorded at the end of appendix B. While still straightforward to evaluate the Mellin transform numerically, the result cannot be written in terms of conventional harmonic sums and it is an interesting open problem to characterize this new class of sums.
Finally, we have compared our result for m = 0 with the recent works [15,16], which exploited, respectively, integrability of the theory and high-loop data in the collinear limit. After converting to our basis, we found perfect agreement with both references (showing in particular that they agree with each other). The coordinate space kernel (5.6), its corresponding eigenvalue for m > 0, and the nonlinear terms in eq. (4.34c), are new predictions.

Collinear singularities and resummation
The eigenvalue is plotted for m = 0 and m = 1 in figures 5-7. It is apparent that, especially near the peak for m = 0, the perturbative series suffers from slow convergence. This was JHEP02(2018)058 Figure 5. The BFKL eigenvalue for m = 0 along the real ν axis at various orders for λ = g 2 YM N c = 6. Convergence near the maximum is visibly slower than away from it. The "resummation of leading-order" is defined below eq. (5.16).
observed already at two loops and explained in terms of nearby singularities in the complex plane at iν = ±1 [12].
In short, these singularities are related to the collinear limit of BFKL, where the scaling dimension ∆ = 2+iν = 3 of the exchanged state coincides with that of twist-two operators: ∆ = 2 + j + γ(j) with j close to 1, e.g. the operators entering the DGLAP equation. As is common for two-level quantum systems, this crossing of two energy levels [18] gets resolved as depicted in figure 6: 16π 2 , one branch choice gives the near-horizontal BFKL trajectory while the other gives the 45 • twist-two (DGLAP) trajectory. (The square root formula follows easily by solving ∆(j) ≈ j + 2 + 8g 2 j−1 for the j, within the overlapping regime of validity of BFKL and DGLAP g 2 |j − 1| 1 where the anomalous dimension γ(j) can be approximated by its leading pole.) It was shown that, expanding the square root to order g 4 , reduces by half the magnitude of the two-loop corrections to the intercept j(0, 0) (if one also includes the complex conjugate singularity at iν = −1) [12]. The "LO resummation" curve in figure 5, called "scheme 2" in ref. [12], thus shows the LO trajectory plus eq. (5.16) minus its O(g 2 ) expansion. (It would be useful to develop a NLO resummation and we leave it as an open problem for the future.) The formula (5.16), expanded to three loops, turns out to not predict very well the three-loop correction to the intercept j(0, 0) ≈ 1 + 11.09g 2 − 84.08g 4 − 2543.05g 6 + O(g 8 ).
In fact it gets even the sign wrong. By looking at the singular terms in F close to the pole we can try to understand why: where δ = 1 − iν. Comparing with eq. (5.16), we find that the leading pole 1024g 6 /δ 5 is exactly as predicted (as it had to). Setting δ = 1, the subleading poles however also give a numerically large contribution to the intercept 2F , so truncating to the leading pole does not give a good approximation to the intercept. However, summing up all the singular terms in eq. (5.17), one finds that about 80% of the three-loop correction to the intercept is reproduced. A heuristic explanation is that the contributions from the next singularities, at iν = ±3, are suppressed by their distance.
Interestingly, all polar terms at L-loops can be obtained from the L-loop DGLAP equation. (See for example [59,60].) From the higher-loop DGLAP equation one can get JHEP02(2018)058 nonsingular terms in the expansion (5.17), see for example eq. (21) of [16]. We have verified that our result (5.15)-(C.3) agrees with all these constraints. 7 We conclude that the physical picture of [12], that large corrections to the intercept originate from the iν = ±1 collinear singularities, is consistent with the three-loop trajectory we obtained, although the full polar part, predicted by DGLAP (as opposed to just the leading pole), must be retained. In general it would be very interesting to find a way to make full use of the DGLAP information at a given loop order Finally, we comment on the Mellin transform of the level-crossing formula (5.16) back to coordinate space. The transform produces a Bessel function: The right-hand side has appeared in coordinate space and momentum space resummations [13,14], so it is nice to see how it arises form the familiar two-level crossing formula (5.16).

Discussion and conclusion
In this paper we have computed, for the first time, the evolution equation which resums large rapidity logarithms in forward scattering to three loops in a gauge theory. Our main results are the full nonlinear equation (4.34) in planar supersymmetric Yang-Mills theory, its linearization (5.6), characterizing the BFKL Pomeron in impact parameter space, as well as its eigenvalue, the Pomeron Regge trajectory, described in appendix C. This result is a first step toward the analogous QCD result, and by itself can already be used to assess the convergence of perturbation theory and its proposed resummations, and shed light on nonlinear saturation effects at finite coupling. This computation was made possible thanks to a recently established correspondence with the resummation of large so-called non-global logarithms, which occur when soft radiation is excluded from a fixed angular region. This correspondence is helpful because it makes available a body of knowledge on the factorization of infrared and collinear divergences, and at a conceptual level it defines in a clear way the evolution equation to all loop orders. This allowed us to derive a systematic subtraction method for nested subdivergences, embodied in eq. (4.21), such that all energy integrals at fixed angle become convergent. We then dealt with collinear subdivergences and real-virtual cancellations in a second step, by multiplying and dividing by the corrections to the single soft current as in eq. (4.14).
Therefore, although we set up our calculation in dimensional regularization and some divergent intermediate objects appeared, we find that in the end the evolution equation depends only on the → 0 limit of physical scheme-independent quantities like the the Bern-Dixon-Smirnov remainder (4.16)! This opens a new possibility to relate an object with the topology of the cylinder, the BFKL Pomeron, to the integrable system appearing in planar scattering amplitudes [61]; graphically speaking, this cuts the cylinder into two half-pipes.

JHEP02(2018)058
As a highly nontrivial test, of both our calculation and of the integrability approach, we have compared our extracted Pomeron trajectory (5.15) with the recent predictions for m = 0 in [15,16], and found perfect agreement! We have also found perfect agreement, in the collinear limit, with the prediction from anomalous dimensions of twist-two operators. The trajectory for other transverse angular momenta m, and nonlinear interactions, are new predictions which it would be very interesting to check within the integrability approach.
It is important to clarify the 1/N c counting in which our result is valid. The projectile is assumed to be made of a finite ∼ N 0 c number of Wilson lines, but whose expectation values across the target can be finite, 1 Nc Tr[U 1 U † 2 ] ∼ 1. This asymmetric setup, motivated for example in proton-nucleus collisions, is the same as that for which the Balitsky-Kovchegov equation is strictly derived. In the context of AdS/CFT this counting would apply to e.g. a light probe of a black hole. This is also a well-defined setup and in fact it would be interesting to work out the nonlinear terms in the Balitsky-Kovchegov equation at strong coupling λ, including perhaps 1/ √ λ stringy effects. The linear terms, which govern correlators of light operators with large but not-to-large energies (before the onset of saturation, such that 1 − U ∼ s j 0 −1 /N 2 c 1) have already been identified with graviton exchange [18]. There are several directions in which this work could be extended. One is to go beyond the planar limit at weak coupling, where the two-loop corrections have recently become available [23,62,63]. Interesting new physical effects appear at three loops in the nonplanar sector, for example the 4→2 reggeon transition which "closes the Pomeron loop" and restores the symmetry between the target and projectile would first be seen there (see for example [19]). Through the KLN theorem, the three-loop evolution could also independently predict from real corrections, and thus test, the recent result for three-loop soft anomalous dimension [64]. Another direction is towards QCD: technically, our setup gives direct access to the evolution equation for non-global logarithms, which in QCD will differ from rapidity evolution by terms proportional to the β-function. These could thus be calculated subsequently by calculating matter loops on both sides of the correspondence.

Acknowledgments
SCH would like to thank Kolya Gromov and Vitaly Velizhanin for discussions regarding the integrability prediction for the trajectory. SCH's work was supported partly by the Danish FNU Grant No. 126152 and by the Danish National Research Foundation (DNRF91). MH is supported by the Villum Foundation Grant No. YIP/VKR022599.

A One-loop correction to the squared double soft current
Here we record the interference of the tree and one-loop double soft current, defined in eq. (4.8), obtained from the soft limit of the six point amplitude as explained in the text. with X = −2 c Γ 2 (Q 2 [1 00 2] /µ 2 ) − + 2π 2 + O( ), and the parity operation P : {1, 0} ↔ {2, 0 }. Here all analytic continuations have been performed, so the logarithms are all real for timelike (positive) invariants, as is the case for our application. The infrared divergences are contained in the factor X but in practice all we will need is the fully renormalized form factor, defined in eq. (4.14), which is finite and obtained by a simple substitution:

B Doing transverse integrals efficiently
Two-dimensional integrals can be done extremely efficiently with the differential equation method. Here we elaborate on our implementation, emphasizing the simplifications related to the fact that all the integrals are absolutely convergent and done directly in 2 dimensions. We first illustrate the method on the integral which occurs at two-loops when obtaining the translation-invariant kernel H(y) (5.10).
The idea is to differentiate with respect to y and add a total derivative with respect to z to remove derivatives of rational factors. Indeed, using the relevant identity: one readily gets that We "win" because the derivatives commutes with the rational factor and hits the logarithm, producing a simpler integral.
An important subtlety is that the left-hand-side of eq. (B.2) is singular and so the equation is only strictly valid for generic z. There are additional contact terms given by the "holomorphic anomaly" This can be understood from the two-dimensional Poisson equation ∂ z ∂z log(zz) = πδ 2 (z). These terms would be absent in dimensional regularization but appear because we insist JHEP02(2018)058 to work with = 0 (see [54] for four-dimensional examples). In the example (B.3), the contact terms are at z = 1 and z = y but the logarithm turns out to vanish on both, so these can be dropped. Evaluating the derivative then gives simply d dy g 1 (y,ȳ) = 1 y d 2 z π (1 + y) z(1 + y − z) . (B.5) To finish, one can repeat the same procedure, inserting a variant of eq. (B.2) to differentiate the integral with respect to y (and/orȳ). Now only the contact term contributes and a general result obtained this way is = log |a − d| 2 |b − c| 2 |a − c| 2 |b − d| 2 , (B.6) which gives (using the vanishing at y = −1 to fix the integration constant) d dy g 1 (y,ȳ) = 1 y × 2 log(yȳ) −→ g 1 (y,ȳ) = log 2 (yȳ) .

(B.7)
This result can be easily confirmed by numerical integration. A critical point to emphasize is that the factor 1/y in eq. (B.5) had to be pulled out in front of the integral before taking the second derivative. If the derivative were allowed to act on that factor, one would gain nothing from it. Only properly normalized integrals simplify upon taking derivatives.
A simple criterion to identify properly normalized integrals is that all the Poincaré residues of their rational factors should be constant (these are often called leading singularities). These are simply the double residue with respect to z followed byz, of the rational factors in the integrand, with z andz treated as independent complex variables. This property is easily verified in eqs. (B.1) and (B.6). Its significance is that it ensures that derivatives of the rational factors have vanishing Poincaré residues, which is needed for them to be total derivatives which simplify upon integration by parts as in eq. (B.2). See for instance refs. [53,54] for other applications of this criterion.
The procedure to decompose an integral into properly normalized ones is essentially partial fractions. When the denominators do not couple z andz, it is in fact literally partial fraction in these two variables, one after the other. But the integrals we need also contain in the denominator an irreducible quadratic form Q(z,z), which is harder to partial-fraction out. We illustrate this with the other integral appearing in H (2) (y), coming from the second term in the kernel (3.12): log |1 − z| 2 |y − z| 2 |1 + y − z| 2 |z| 2 , Q = |1−z| 2 |y −z| 2 −|1+y −z| 2 |z| 2 .
Despite appearances, Q is a quadratic form in z,z. The leading singularities of the rational factor can be computed and found to be linear combinations of 1/[(1 − y)(1 +ȳ)] and 1/[(1 + y)(1 −ȳ)], so the decomposition into properly normalized integrals will require two terms. To illustrate the result of the partial-fraction method to be detailed shortly, one indeed finds that the integral can be rewritten exactly as: be easily fixed from the limit x → 1. Its symbol turns out to be made of the five letters x, 1 − x,x, 1 −x and 1 − v = x +x − xx. At transcendental weight 3, we found only three nontrivial single-valued functions with such symbol, in terms of which the three-loop linearized kernel K (3)lin (x) in (5.6) is compactly written: Although not manifest from these formulas, these functions have no branch cut on the complex plane wherex = x * . This can be confirmed by series-expanding around singular points such as x =x = 0 and x =x = 1, where to all orders one finds only singlevalued logarithms of log(xx) or log(1 − x)(1 −x) but never log(x) nor log(1 − x) separately. Furthermore, there are no singularities along 1 − v = 0 (which traces a unit circle with center at x = 1).
We note that the same five letters are also singularities of the two-loop kernel, so it is natural to conjecture that no other letters appear in K (L)lin (x) to any order in perturbation theory in planar N = 4 SYM. Its translation-invariant projection H (L) (y), defined by the integration (5.10), can then be obtained by applying the algorithm detailed above, which implies that at most the ten letters d log y,ȳ, 1±y, 1±ȳ, y +ȳ, 1 + yȳ, This is regular and in fact vanishes at ν = 0, in accordance with the all-order result (5.8).
Other values of m can be evaluated numerically using the attached Mathematica notebook.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.