Two-parton scattering amplitudes in the Regge limit to high loop orders

Abstract
We study two-to-two parton scattering amplitudes in the high-energy limit of perturbative QCD by iteratively solving the BFKL equation. This allows us to predict the imaginary part of the amplitude to leading-logarithmic order for arbitrary t-channel colour exchange. The corrections we compute correspond to ladder diagrams with any number of rungs formed between two Reggeized gluons. Our approach exploits a separation of the two-Reggeon wavefunction, performed directly in momentum space, between a soft region and a generic (hard) region. The former component of the wavefunction leads to infrared divergences in the amplitude and is therefore computed in dimensional regularization; the latter is computed directly in two transverse dimensions and is expressed in terms of single-valued harmonic polylogarithms of uniform weight. By combining the two we determine exactly both infrared-divergent and finite contributions to the two-to-two scattering amplitude order-by-order in perturbation theory. We study the result numerically to 13 loops and find that finite corrections to the amplitude have a finite radius of convergence which depends on the colour representation of the t-channel exchange.


Introduction
The study of QCD scattering in the Regge limit has been an active area of research for over half a century, e.g. [1][2][3][4][5][6][7]. While the general problem of high-energy scattering is non-perturbative, in the regime where the exchanged momentum −t is high enough, i.e. s −t Λ 2 QCD (see figure 1), perturbation theory offers systematic tools to analyse this limit. Central to this is the Balitsky-Fadin-Kuraev-Lipatov (BFKL) evolution equation [1,2], which provides a systematic theoretical framework to resum high-energy (or rapidity) logarithms, ln(s/(−t)), to all orders in perturbation theory. This approach was used extensively to study a range of physical phenomena including the small-x behaviour of deep-inelastic structure functions and parton densities, and jet production with large rapidity gaps. Furthermore, non-linear generalisations of BFKL, known as the Balitsky-JIMWLK equation [8][9][10][11][12][13], are today a main tool in the theoretical description of dense states of nuclear matter, notably in the context of heavy-ion collisions.
While many applications of rapidity evolution equations to phenomenology require the scattering particles to be colour-singlet objects, in the present paper we are concerned with the more theoretical problem of understanding partonic scattering amplitudes in the high-energy limit, similarly to refs. [14][15][16][17][18][19][20][21][22][23][24][25]. This is part of a more general programme of understanding the structure of gauge-theory amplitudes and the underlying physical and mathematical principles governing this structure. The basic observation is that gauge dynamics drastically simplifies in the high-energy limit, which renders the amplitudes computable to all orders in perturbation theory, to a given logarithmic accuracy.
The present paper continues our recent study [23][24][25] of 2 → 2 partonic amplitudes (qq → qq, gg → gg, qg → qg) in QCD and related gauge theories. A key ingredient in these studies is provided once again by rapidity evolution equations, BFKL and its generalisations, which are used to compute high-energy logarithms in these amplitudes order-by-order in perturbation theory.
Scattering amplitudes of quarks and gluons are dominated at high energies by the tchannel exchange (figure 1) of effective degrees of freedom called Reggeized gluons. 2 → 2 amplitudes are conveniently decomposed into odd and even signature characterising their symmetry properties under s ↔ u interchange, or crossing symmetry: where odd (even) amplitudes M (−) (M (+) ) are governed by the exchange of an odd (even) number of Reggeized gluons. Furthermore, as shown in ref. [24], these have respectively real and imaginary coefficients, when expressed in terms of the natural signature-even combination of logarithms, The real part of the amplitude, M (−) , is governed, at leading logarithmic (LL) accuracy, by the exchange of a single Reggeized gluon in the t channel. To this accuracy, high-energy logarithms admit a simple exponentiation pattern, namely where the exponent is the gluon Regge trajectory (corresponding to a Regge pole in the complex angular momentum plane), α g (t) = αs π C A α (1) g (t) + O(α 2 s ), whose leading order coefficient α (1) g (t) is infrared singular, α (1) g (t) ∼ 1 2 in dimensional regularization with d = 4−2 (see eq. (2.3) below). Infrared singularities are well-known to exponentiate, independently of the high-energy limit. Importantly, however, eq. (1.3) illustrates the fact that the exponentiation high-energy logarithms must be compatible with that of infrared singularities, which is a nontrivial constraint on both. This observation and its extension to higher logarithmic accuracy underpins a long line of investigation in refs. [14][15][16][17][18][19][20][21][22][23][24][25].
The key property of the Reggeized gluon being signature-odd greatly constrains the structure of higher-order corrections. For the real part of amplitude, the simple exponentiation pattern generated by a single Reggeized gluon is preserved at the next-to-leading logarithmic (NLL) accuracy, except that it requires O(α 2 s ) corrections to the trajectory and also the introduction of (s-independent) impact factors. This simple picture only breaks down when three Reggeized gluons can be exchanged, which first occurs at NNLL accuracy and leads to Regge cuts. This contribution was computed in ref. [24] through three-loops, by constructing an iterative solution of the non-linear Balitsky-JIMWLK equation which tested the mixing between one and three Reggeized gluons.
In this paper we focus on the imaginary part of the amplitude, M (+) , extending our work [25]. Here the leading tower of logarithms, in which we are interested, is generated by the exchange of two Reggeized gluons, starting with a non-logarithmic term at one loop: Here we suppressed subleading terms in as well as multiloop corrections, which take the form α s L −1 at loops; because the power of the energy logarithm L is one less than that of the coupling, these are formally next-to-leading logarithms (NLL). In eq. (1.4) one may observe another salient feature of this tower of corrections, namely the colour structure, which is even under s ↔ u interchange (M tree is odd, and so is the operator T 2 s−u acting on it). The first term in the square brackets in (1.4) is the exact result in the planar limit; we will be interested in the full series of corrections α s L −1 , which are all subleading in the large N c limit (see the definitions of colour operators in eq. (2.8) below).  they are particularly important, and clearly at high energies, where α s L ∼ O(1), one should aim at an all-order calculation. These corrections, however, do not feature a simple exponentiation pattern as in eq. (1.3); they give rise to a Regge cut rather than a pole. We shall study these corrections using an iterative solution of the BFKL equation, continuing the work of ref. [23][24][25]. In [23] higher-order terms in eq. (1.4) were computed through four loops -the first order where finite contributions appear (see eqs. (28)(29) in [24]). Subsequently, in ref. [25] infrared-singular contributions were computed in dimensional regularization to all orders. The purpose of the present paper is to extend the calculation to finite contributions, and in particular, to obtain the infrared-renormalized amplitude, or hard function, which we expect (together with the soft anomalous dimension) to control any infrared-safe cross section. We are interested in the exact perturbative solution of the BFKL equation for any colour exchange, that is, not restricted to the planar limit. While the BFKL Hamiltonian was famously diagonalized by its authors in the case of color-singlet exchange, the solution is not known in the general case. Adding to the complexity is the fact that amplitudes are infrared singular, forcing us to work in dimensional regularization. While it is not known how to diagonalise the BFKL Hamiltonian in these circumstances, we are able to solve the problem by using two complementary approaches, the first by taking the soft approximation while maintaining dimensional regularization, and the second by considering general (hard) kinematics in strictly two transverse dimensions. Let us briefly describe each of these approaches.
The first approach is a computation of the wavefunction describing the emission of two Reggeons at ( − 1) loops, and the corresponding -loop 2 → 2 amplitude, in the soft approximation, where one of the two Reggeized gluons carries transverse momentum k 2 which is significantly smaller than the total momentum transfer by the pair, −t = p 2 , i.e. the limit characterized by a double hierarchy of scales k 2 p 2 s. This is the limit used in ref. [25] to determine all infrared-singular contributions the amplitude. This was achieved using the simple observation that the wavefunction is itself finite to all orders in perturbation theory and that BFKL evolution closes within this approximation. All the singularities of the amplitude at any given loop order are in turn produced in the final integration over the wavefunction (corresponding to the transition from the top to the bottom row in figure 2). In the present paper, building upon the computation of the wavefunction in [25] we introduce a symmetrized solution accounting simultaneously for the two soft limits, k 2 p 2 and (p−k) 2 p 2 , which amounts to an elegant separation between soft and hard contributions to the wavefunction and amplitude. Within this approximation we are able to write down a resummed analytic expression for the amplitude, including its finite contributions.
The second approach, which we develop in the present paper, is based on starting with the BFKL equation in exactly two dimensions. Without making any further approximation, we set up an iterative solution of the equation by identifying differential operators that commute with (parts of) the Hamiltonian up to a computable set of contact terms. Evolution induced by the Hamiltonian then becomes trivial within a class of iterated integrals dictated by the nature of the problem, these are the Single-Valued Harmonic Polylogarithms (SVHPLs), first systematically classified by Francis Brown in ref. [26] and then studied and applied in the context of motivic periods [27] and Feynman integrals [28,29]. The relevance of this class of functions for gauge-theory amplitudes within the Regge limit [30][31][32][33][34][35] (and beyond [36,37]) has been recognised in recent years, and it is important also in our current problem: the hard wavefunction, defined in strictly two dimensions, is fully expressible in terms of SVHPLs, and the corresponding contribution to the amplitude can in turn be written in terms of Single-Valued Multiple Zeta Values (SVMZVs). For the ladder graphs relevant here, each additional loop increases the transcendental weight by one unit. The resulting uniform-weight expressions in terms of single-valued functions are significantly simpler as compared to the corresponding ones in terms of ordinary polylogarithms and zeta values. For the final integration over the wavefunction we develop two independent approaches, one relying on analytic continuation and integration over the discontinuities of the wavefunction away from the region were they are single-valued, and the other relying instead on a modified application of the evolution algorithm itself. The two yield identical results. By combining the hard contribution to the amplitude with the dimensional-regularized soft contribution we compute the full amplitude, in principle to any order, and in practice to thirteen loops.
The structure of the paper is as follows. In section 2 we present the BFKL equation in dimensional regularisation, bring it to a form suitable for iterative solution and review the relation between the off-shell wavefunction and the two-to-two scattering amplitude. We also show how an iterative solution can be obtained for the first few orders directly in dimensional regularization without resorting to any approximation, and explain why this approach does not practically extend to higher orders. In this context we compute the amplitude numerically through five loops, providing a valuable check for our subsequent calculations. Next, in section 3 we review the soft approximation developed in [25] and explain how infrared factorization, combined with the finiteness of the wavefunction, facilitate a systematic separation of the latter into 'soft' and 'hard' components, such that eventually, finite corrections to the infrared-renormalized scattering amplitude can be determined in full. To this end we introduce a symmetrized version of the soft wavefunction, which captures both soft limits, and then derive an analytic expression for the amplitude as a function of α s L, which resums both infrared-divergent and finite contributions to all loops, within the soft approximation. In section 4 we turn to discuss the wavefunction in general (hard) kinematics. Working directly in two dimensions we introduce the relevant kinematic variables, analyse the action of the BFKL Hamiltonian and demonstrate that evolution generated by this Hamiltonian translates into an algorithmic procedure in the space of SVHPLs. Having determined the wavefunction order by order, we turn in section 5 to compute the corresponding two-to-two scattering amplitude. In section 6 we perform a numerical study of the resulting wavefunctions and amplitudes, and address the convergence of the perturbative expansion. Finally, in section 7 we make some concluding comments and present an outlook for future investigation.

BFKL equation in dimensional regularisation and the → amplitude
In the high-energy limit, scattering amplitudes are conveniently described in terms of Wilson lines, which dress the external partons. The evaluation of vacuum expectation values of Wilson lines stretching from minus to plus infinity leads to rapidity divergences, which needs to be renormalised. As a consequence, the renormalised amplitude obeys a rapidity evolution equation, which can be shown to correspond to the Balitsky-JIMWLK equation. In this paper we are interested to study the two-Reggeon exchange contribution to two-parton scattering amplitudes, for which the evolution equation reduces to the BFKL equation [23,24]. The scattering amplitude can be determined formally to any order in perturbation theory as an iterative solution of the dimensionally-regularised BFKL equation. This procedure was described in [25], to which we refer for further details. In this section we review the definitions necessary to set up the calculation.
In the following we consider the two-reggeon exchange contribution to 2 → 2 scattering amplitudes. We can single out this contribution by introducing a reduced amplitude, in which the one-Reggeon exchange has been removed: where L is the signature-even high-energy logarithm defined in eq. (1.2), T 2 t represents the total colour charge exchanged in the t channel (see eq. (2.8) below) and i, j are the species indices defining the two-parton scattering; in what follows we will drop these indices, unless explicitly needed. Finally, the function is the gluon Regge trajectory introduced already in eq. (1.3), where the leading-order coefficient in dimensional regularization with d = 4 − 2 is given by belongs to a class of bubble integrals which will be defined below.
The two-Reggeon cut contributes only to the even amplitude defined in eq. (1.1), thus we focus only on this component in the following. As discussed in [25], the reduced amplitude takes the form of an integral over the two-Reggeon wavefunction Ω(p, k), as follows: where p 2 = −t. In eq. (2.5) the integration measure is and M (tree) ij→ij represent the tree amplitude, given by where λ i for i = 1 through 4 are helicity indices. The colour operator T 2 s−u in eq. (2.5) acts on M (tree) ij→ij and it is defined in terms of the usual basis of quadratic Casimirs corresponding to colour flow through the three channels [22,38]: where T i is the colour-charge operator [39] associated with parton i. The BFKL equation [1,2] for the wavefunction Ω(p, k) in eq. (2.5) takes the form where L is the high-energy logarithm (1.2) and where the Hamiltonian takes the form [25] where two independent colour factors come along with two different operations: 11b) The function f (p, k, k ) in eq. (2.11a) represents the BFKL evolution kernel and J(p, k) in eq. (2.11b) is defined by While it is unknown how to diagonalise this d-dimensional Hamiltonian, we may invoke a perturbative solution [23,25] by expanding the wavefunction in the strong coupling constant: where we set the renormalisation scale equal to the momentum transfer, µ 2 = −t = p 2 . Substituting the expanded form of the wavefunction in (2.14) into the BFKL evolution equation (2.9) one deduces that whereĤ is the BFKL hamiltonian of eq. (2.10), that is, the wavefunction at any given order is found by repeated application of the BFKL Hamiltonian, where the initial condition in our normalization is simply Next, let us consider the on-shell 2 → 2 amplitude. Substituting the expanded wavefunction (2.14) into (2.5) we readily obtain the following expansion Namely, integrating over the ( − 1)-th order contribution to the wavefunction yields the -th order contribution to the amplitude.
A graphical illustration of eq. (2.18) is provided in figure 3. As discussed in the introduction, because of BFKL evolution, the amplitude at NLL accuracy can be represented as a ladder. At order it is obtained by closing the ladder and integrating the wavefunction of order ( − 1) over the resulting loop momentum, according to eq. (2.18). The wavefunction Ω ( −1) (p, k) in turn is obtained by applying once the leading-order BFKL evolution kernel to the wavefunction of order ( − 2). Graphically, this operation corresponds to adding one rung to the ladder.
. Graphical representation of the amplitude at NLL accuracy, as obtained through BFKL evolution. The addition of one rung corresponds to applying once the leading-order BFKL evolution on the wavefunction of order ( − 2). This gives the wavefunction at order ( − 1), according to eq. (2.10). Closing the ladder and integrating over the resulting loop momentum gives the reduced amplitude, according to eq. (2.18).
Inspecting eqs. (2.11a) and (2.11b) we see that the BFKL evolution consists of an integration and a multiplication part. The effect of evolution is thus expressed formally in a compact form by introducing a class of functions where Ω ∅ (p, k) ≡ 1, and w indicates a word made of indices "i" or "m", which stand for integration and multiplication, respectively, according to the action of the two Hamiltonian operators in eq. (2.11a) and (2.11b), respectively. In this notation the first four orders of the wavefunction read, for instance, ) ) Symmetries play an important role in determining the general structure of the wavefunction, and from a practical perspective they can be useful to reduce the number of integrals that need to be evaluated at each loop order. The wavefunction is symmetric under swapping the two t-channel Reggeons, which can be understood from the graphical representation of the BFKL evolution in figure 3. This implies which can be easily verified by showing that the functions f (p, k, k ) in (2.12), J(p, k) in (2.13) and Ω (0) (p, k) in (2.16) obey the same symmetry. This symmetry property will become handy in section 3, making it possible to capture simultaneously both soft limits, k 2 → 0 and (p − k) 2 → 0. This, in turn, will be important for implementing a systematic separation between the soft and hard regimes, without needing an extra regulator. Despite the simplifications allowed by symmetries, though, the evaluation of the wavefunction in 2 − 2 transverse dimensions without additional simplifications becomes quickly infeasible. For instance, already the wavefunctions with one or two integrations (one or two occurrences of the index "i") involve integrals of the type which are represented respectively in figure 4 (a) and (b). Such integrals evaluates to Appell, and more in general Lauricella functions in dimensional regularisation. Given the lack of a systematic classification of these functions in terms of iterated integrals, the evaluation of the wavefunction beyond the third order is not practical.  . Three-mass triangle integrals with massless propagators, which appear in the calculation of the wavefunction at two and three loops. These integrals contribute to the amplitude only starting respectively at four and five loops, due to symmetry constraints, as discussed in the main text. The bubble integral on one of the edges of the triangle clarifies the origin of the propagator which is raised to power in eq. (2.25).
The amplitude at order is obtained upon integrating the wavefunction of order − 1, as indicated in eq. (2.18). As in case of the wavefunction, symmetries turn out to be important for a simplification of the calculation and interpretation of the result. While the two Reggeons in the wavefunction can be defined to originate from either the projectile or target Wilson line -which gives the corresponding ladder graphs a sense of direction -this is no longer true at the level of the amplitude. Physically the two cases become indistinguishable, and we refer to this as the target-projectile symmetry. In general, this implies the relation [25] [Dk] (2.26) Furthermore, in the notation of eqs. (2.19a) and (2.19b) reversal of the rungs directly translates to the reversal of the indices of the wavefunction. The target-projectile symmetry thus guarantees the equality (2.27) The symmetries discussed above can reduce the number of functions to be computed significantly, and make the calculation of the amplitude trivial up to three loops, since it can be shown that the integration of the wavefunction involves only bubble integrals. Furthermore, the calculation of the amplitude at four loops in dimensional regularisation is still feasible, as it involve bubble integrals and a single more involved kite-like integral, represented in figure 5 (a). Up to four loops one obtains [25] M (+,1) A thorough discussion of the target-projectile symmetry, and its effect on the colour structure of the amplitude has been given in [25], to which we refer for further details. In this paper we are interested to evaluate the amplitude, including finite terms, to higher orders in the perturbative expansion. Despite the symmetries discussed above, however, beyond four loops the iterated integrals appearing are all but easy with current methods. A simple and fast way to extend the study in ref. [23,25] to higher loops is provided by numerical integration methods. In particular, we find sector decomposition as implemented in pySecDec/SecDec [40,41] to be suited to calculate the nested integrals that enter the five-loop amplitude. Provided a high numerical accuracy it is straightforward to extract from the results the rational coefficients of the zeta numbers appearing at this loop order. This procedure relies on the observed homogeneous transcendental weight property of the -loop amplitude: Assigning o( ) = −1, o(π) = 1 and o(ζ n ) = n one sees that the terms of the -loop amplitude are uniformly of weight o(M (+, ) NLL ) = . We can hence deduce which zeta numbers (or powers of π) may appear at any given order in .
Another observation facilitates this procedure at five loops; after dividing the -loop amplitude by B 0 (2.4) there are no occurrences of ζ 2 = π 2 /6 up to four loops, see e.g. the O( ) terms of eq. (2.31). If we assume this absence of ζ 2 to be an actual property of the amplitude, the finite terms of the five-loop amplitude can only be proportional to one transcendental number, ζ 5 , whereas ζ 3 ζ 2 is excluded. At this point this approach may seem rather conjectural. However, over the course of the next two sections we develop methods that prove this assumption, and we shall briefly return to it at the end of section 5.3.
To obtain the five-loop amplitudeM (+,5) NLL we integrate the four-loop wavefunction Ω (4) (p, k) of (2.23) according to eq. (2.18). In doing so one is faced with a plethora of multi-loop integrals. Many of them correspond to bubble graphs and can be easily evaluated analytically. Others vanish because of the symmetries discussed above. The remaining integrals can be computed numerically using pySecDec. One of the more difficult examples is shown in figure 5. In the depicted case one can integrate out the two internal bubbles and is left with a three-loop integral with two of the propagators raised to non-integer powers: (2.32) After combining all contributions (and reconstructing the zeta numbers in case of the (a) (b) p p Figure 5. Example of a four-and five-loop integrals that enters the calculation of the four-and five-loop amplitude respectively. The two bubbles may be integrated out, turning it into a two-and three-loop integral with two propagators raised to non-integer powers, cf. eq. (2.32).
numerical results) we find This result will serve as a consistency check for our computation below.

The soft approximation
In section 2 we have shown how the two-Reggeon contribution to the two-parton scattering amplitude is conveniently described in terms of the reduced amplitudeM. The latter is defined in eq. (2.1) by (multiplicatively) removing the single-Reggeon effect from the full amplitude M. This allowed us to use BFKL evolution to express the two-Reggeon contribution toM in terms of iterated integrals. Beyond four loops these integrals become difficult to evaluate exactly in d = 4 − 2 dimensions, but as we are going to show now, this is also not necessary. Ultimately we are interested in extracting physical information about the scattering process, and dimensional regularization is used in the present context for the sole purpose of regularizing long-distance singularities 1 . Here infrared factorization come into play: the long-distance singularities of M can be factorized, M = ZH, where the "infrared renormalization" factor Z captures all divergences (which famously exponentiate in terms of the soft anomalous dimension, see e.g. [17,36,39,[42][43][44][45][46][47][48][49]) while the infrared-renormalized amplitude H -sometimes referred to as the "hard function" -is finite, and can be evaluated in four space-time dimensions (or equivalently, two transverse dimensions). To understand this from a physical perspective recall that physical quantities such as cross sections are finite: starting from the infrared-singular amplitude M, their calculation inevitably incorporates a mechanism of cancellation of the singularities involving soft real-gluon emission. Once this was implemented, the finite, physical result can only depend on four-dimensional quantities, namely the soft anomalous dimension and the infrared-renormalized amplitude H.
In Ref. [25] we have shown that the soft anomalous dimension associated with the signature-even amplitude, or indeed the relevant infrared renormalization factor Z, can be computed to all orders by evaluating the reduced amplitudeM to O( −1 ). Similarly, we are going to show now (section 3.1) that the infrared-renormalized amplitude H (in four dimensions) can be completely determined from the reduced amplitudeM, evaluated at the same accuracy, i.e. to O( 0 ). This, along with the fact that the corresponding wavefunction Ω is finite, greatly simplifies the task of performing BFKL evolution to high loop orders, because it allows us to follow an "expansion by region" approach: in section 3.2 we split the wavefunction into soft and hard components, each of which is rendered computable using different considerations. The soft wavefunction -giving rise to all the singularities in the amplitude -can be computed analytically in dimensional regularization owing to the drastic simplification of BFKL evolution in this limit, while the hard wavefunction is only required in strictly two transverse dimensions, where BFKL evolution again simplifies (see section 4). These two wavefunction components will subsequently serve to compute the corresponding soft and hard contributions to the reduced amplitudeM to the required order, O( 0 ). In section 3.3 we review the main results of Ref. [25] regarding the all-order computation of the wavefunction within the soft approximation. We also introduce there a symmetrized soft wavefunction which captures both soft limits. This, in turn, is used in section 3.4 to compute the corresponding O( 0 ) contributions to the reduced amplitude.
Finally, in section 3.5 we make use of the results of sections 3.1 and 3.4 to evaluate the O( 0 ) soft contributions to the infrared-renormalized amplitude H.

Infrared factorisation in the high-energy limit
According to the infrared factorisation theorem (see e.g. [17,36,39,[42][43][44][45][46][47][48][49]), infrared singularities of an amplitude M are multiplicatively renormalised by a factor Z, such that the infrared-renormalized amplitude H is finite as → 0. We use a minimal subtraction scheme, where the renormalisation factor Z consists of pure poles. It is then given explicitly as the path-ordered exponential of the soft anomalous dimension: where, to the accuracy needed in this paper, we can restrict to tree-level running coupling: α s (λ) = α s (p) p 2 /λ 2 . Given that Z was determined in Ref. [25] to NLL accuracy in the high-energy logarithm, our goal here is to determine the infrared-renormalized amplitude H to the same accuracy. Thus we need to specialise eq. (3.1) to the high-energy limit. Recalling that in this limit the amplitude splits naturally into even and odd components under the s ↔ u signature symmetry, we may focus directly on the even component (the odd component was analysed already in [24]): Our final goal is to determine H where we defined x ≡ αs π L. The factor Z (−) NLL was determined to all orders in perturbation theory in [25]: comparing eqs. (4.12), (4.14) and (4.17) there we express Z where the function R( ) reads see also eq. (3.16) of [25]. In eq. (3.6) the factor exp[−(xT 2 t )/(2 )] on the l.h.s. is left there, because it corresponds directly to the factor Z −1 (+) LL appearing in eq. (3.4) to the left of Z (−) NLL . Notice also that the −1 in the numerator of the first fraction on the r.h.s. of eq. (3.6) can actually be removed, given that we need to consider only the poles originating from eq. (3.6), and this term contributes only at O( 0 ), given that the second -dependent factor is regular, i.e.
where the factor e − x NLL of (3.6) can be readily substituted as well (this will be done in section 3.5). Eq. (3.9) is an important step because (given that B 0 ( ) − 1 = O( 2 ), eq. (2.4)) it clearly shows that the hard function H

Soft and hard wavefunction and amplitude
Our strategy to compute the finite part of the reduced amplitudeM (+) NLL at higher orders is to separate soft and hard components of the wavefunction and truncate the latter to two transverse dimensions ( = 0), where BFKL evolution is much more tractable (see section 4).
As demonstrated in ref. [25], the soft limit of the wavefunction, where one of the two Reggeons has a small momentum, e.g., k 2 (p − k) 2 p 2 , fully determines all the singular parts in . This was used to obtain the all-order result for the renormalisation factor Z (−) NLL in eq. (3.6). In addition, the soft limit generates some O( 0 ) finite contributions, which must be added to those generated by the complementary hard region, where both k 2 and (p − k) 2 are of order p 2 .
To control O( 0 ) terms a clear separation between the two regions is necessary. We choose to do this at the level of the wavefunction Ω(p, q). Recall that Ω(p, q) is a finite function 2 of [25], i.e. any singularities in the reduced amplitude are generated through the final integration over the wavefunction in (2.5). To proceed we split the wavefunction into two terms: Ω(p, k) = Ω s (p, k) + Ω h (p, k) , (3.10) such that the second term, the hard component, vanishes in soft limits: It then follows from (2.5) that no singularities can be generated upon integrating Ω h (p, k) (i.e. all the singularities inM NLL are generated upon integrating Ω s (p, k)) and hence only the → 0 limit of Ω h contributes to the finite part of the reduced amplitude. Denoting the wavefunction in this limit as Equations (3.13) and (3.14) are central to our approach and will guide our computations in what follows. They show that, to compute the finite part of the reduced amplitude, we must treat the soft wavefunction exactly as a function of , but we are allowed to truncate the hard wavefunction to O( 0 ). Note that in (3.14b) we have already substituted the twodimensional limit of the hard wavefunction, so taking the → 0 limit simply amounts to taking the integration momentum k to be two-dimensional. These finite integrals will be done in section 5. Let us briefly summarise our plan for the reminder of this section. After reviewing the main arguments of [25], our aim in section 3.3 is to present a symmetrized version of the soft wavefunction in dimensional regularization, eq. (3.22), which simultaneously captures the two regions where either of the two Reggeons is soft. We then extract the O( 0 ) terms in the wavefunction and resum them; these will be used in section 5 to determine the two-dimensional hard wavefunction Ω (2d) h (p, k) from the full one according to eq. (3.12). Subsequently in section 3.4 we use the soft wavefunction, computed to all-orders in , to determine the corresponding contributions to the reduced 2 → 2 amplitude. We also present an analytic formula resumming these corrections in eq. (3.36). Finally, in section 3.5 we determine the soft wavefunction contribution to the infrared-renormalized amplitude H (+) NLL using eq. (3.9).

The soft wavefunction
The central property of the wavefunction Ω(p, k) highlighted in [25] and already mentioned above, is the fact that it is finite for → 0, to all orders in perturbation theory. This has far reaching consequences, because it means that all singularities in the amplitude must arise from the last integration in (2.5), and originate from the soft limits k → 0 and k → p of Ω(p, k). One finds that it is particularly easy to calculate the wavefunction in these limits: as it turns out, the soft approximation is closed under BFKL evolution, i.e., starting with Ω (j) (p, k), with k soft, implies that the momentum k in Ω (j−1) (p, k ), which has one rung fewer, can also be taken soft, k → 0, without affecting the result for Ω (j) (p, k). In other words, starting with Ω (j) (p, k) where k is soft is equivalent to considering the entire side rail of the ladder consisting of soft momenta, k , k , . . . → 0. Similarly, starting with k → p implies that all momenta (p−k), (p−k ), (p−k ), . . ., are soft. The symmetry of eq. (2.24), then, implies that Ω(p, k) in the two limits k → 0 and k → p must be the same.
In the soft limit the BFKL hamiltonian becomes [25] is the soft approximation of eq. (2.13). One finds that the wavefunction becomes a polynomial in ξ ≡ (p 2 /k 2 ) , i.e., the soft limit turns BFKL evolution into a one-scale problem. The integrals involved in eq. (3.15) are simple bubble integrals of the type where the integration measure is given in eq. (2.6), and the class of bubble functions B n ( ) is Note that B 0 of (2.4) appearing in the gluon Regge trajectory and in the measure (2.6) corresponds to the special case of (3.18) with n = 0. Using eq. (3.17) one can write the action of the soft Hamiltonian (3.15) on any monomial where we have introduced the notation By making repeated use of eq. (3.19) one finds that the wavefunction at order ( − 1) can be expressed in a closed-form, as follows [25]: .
(3.21) As discussed in [25], this expression can be easily integrated, obtaining an expression for the amplitude which correctly describes its singular part to all orders in perturbation theory.
While eq. (3.21) is perfectly valid in the soft limit, it breaks explicitly the symmetry of eq. (2.24) between the two soft limits. As we will see below, it is advantageous to work with expressions where this symmetry is manifest. In this paper, we thus introduce a different soft wavefunction, obtained by symmetrising eq. (3.21) under k ↔ (p − k): .

(3.22)
This formula simultaneously captures the correct behaviour of Ω(p, k) in both soft limits k → 0 and k → p. It will be used in section 3.4 below to compute the soft contributions to the reduced 2 → 2 amplitude. Before doing that let us have a closer look at the expansion of the soft wavefunction we obtained. We recall [25] that all the negative powers of in (3.22) cancel upon performing the sum over n, leading to a finite wavefunction at any loop order. While positive powers of in (3.22) do play a role in the computation of the amplitude, the leading O( 0 ) have a special role: according to eq. (3.12) it is precisely what must be subtracted from the full two-dimensional wavefunction to obtain the hard wavefunction Ω (2d) h . With this in mind, let us write down explicitly the leading terms in in the first few orders of the soft wavefunction in (3.22): ) ) ) In fact, these terms exponentiate and can be resummed into the following all-order expression using (2.14) for = 0, yielding

Soft contributions to the 2 → 2 amplitude
Next, let us consider the soft contribution to the reduced 2 → 2 scattering amplitudeM. It is straighforward to insert eq. (3.22) into eq. (2.18), perform the last integration and derive the -th order contribution to the amplitude. In particular, given the symmetrised form of eq. (3.22), the last integration can be done with the integration measure [Dk] in eq. (2.6), i.e. avoiding the need to introduce a cut-off as in ref. [25]. After some arrangement we get where the functions B n ( ) andB n ( ) have been defined respectively in eqs. (3.18) and (3.20), and we have introducedB The coefficientsM (+, ) NLL,s in (3.25) are of course polynomial in the colour factors. For illustration, we expand eq. (3.25) to the first few orders in perturbation theory, obtaininĝ NLL,s = iπ NLL,s = iπ NLL,s = iπ where we used the shorthand notation for the colour factors, We note that the expansion coefficients display uniform transcendental weight (where, as usual 1/ has weight 1) and involve exclusively single zeta values (sometimes referred to as ordinary zeta values, namely the values of the Riemann zeta function at integer arguments). We further notice that ζ 2 (or ζ 2 times other zeta values, e.g. ζ 2 ζ 3 at weight 5, etc.) factors do not appear in eqs. (3.27a)-(3.27h) (ζ 2 terms would be present if we were to expand the factor factor B 0 ( )). Higher even zeta numbers do appear, but we will see below that they have a distinct origin as compared to the odd ones.
Given that the expansion coefficientsM (+, ) NLL,s involve just single zeta values, and are moreover of uniform weight, it is interesting to explore the possibility to sum up the series to all orders. Indeed, such summation was achieved for the singular terms in Ref. [25], so let us compare eq. (3.25) above with the result obtained in [25]. There we proved that the singular terms of the reduced amplitude admit a simplified form (3.28) The latter, however, differs from the original soft amplitude obtained from eq. (3.22) starting at O( 0 ) (compare eqs. (3.13) and (3.15) of [25]). A nice feature of eq. (3.28) is that the loop functions in the amplitude at order do not depend on the index , apart from the factor (B 0 ( )) / !, and this allows one to easily resum eq. (3.28) to all orders in perturbation theory, obtaining an expression for the integrated soft amplitudeM NLL,s simpl. : with x = L α s /π (see also eq. (3.18) of [25]). This formula however does not correctly capture the non-singular terms obtained with a cut-off, for which no similar simplification was found. We nevertheless show that an all-order resummation formula can be found for the O( 0 ) corrections to the amplitude defined in our current symmetric scheme, eq. (3.25).
To this end we consider the coefficients defined as the finite part of the difference between those in soft amplitude, eq. (3.25), and in its simplified version, eq. (3.28): After some arrangement the coefficients∆ NLL can be put into the form where we discarded powers of B 0 ( ), which do not affect the finite terms. From eq. (3.31) the coefficients δ ( ) of (3.30) can be determined explicitly in terms of odd ζ numbers and the ratio of colour factors r = such that the sum: We conclude that the series∆ Using now the fact that that the simplified amplitudeM NLL,s simpl. in eq. (3.28) exponentiates independently, see eq. (3.29), we obtain where in the second line we expressed the amplitude in terms of the function R( ) . Writing the reduced amplitude as in the second line of eq. (3.36) makes it easier to extract the infrared-renormalized amplitude from the reduced amplitude, as we will see in section 3.5. Writing explicitly the factor∆ (+) NLL , the reduced amplitude readŝ (3.37) Of course, upon expansion (3.37) yields back the coefficients of (3.25) we listed in (3.27a) through (3.27h). Having at hand a resummed expression we can gain further insight on number-theoretical features of the expansion coefficients in eqs. (3.27a)-(3.27h). We already know based on the derivation above that the∆ NLL component in (3.36) gives rise to odd zeta values only. It then transpires that the sole origin of even ones is the function R( ) in the first term. Further number-theoretical features will be discussed in section 5, once we have computed the hard contribution to the reduced amplitude. The possibility to resum the series for the amplitude to all orders including finite O( 0 ) terms is highly nontrivial, and it is an additional advantage of k ↔ (p − k) symmetric scheme we adopted here for the soft approximation. It will be used below in deriving a resummed expression for the contribution of the soft region to the infrared-renormalized amplitude.

From the reduced amplitude to the infrared-renormalized amplitude
Now that we have determined the soft wavefunction and the corresponding reduced amplitude, we are in a position to consider again the infrared-renormalized amplitude, as defined in eq. (3.9). Following eqs. (3.10) and (3.13) we split the infrared-renormalized amplitude into a soft and a hard component: Then, from eq. (3.9) it follows that where in (3.39b) we neglected positive powers of originating in the expansion of (1−B 0 ( )), using the fact thatM NLL,h is itself finite. Of course such a simplification cannot be applied to (3.39a) where there is an interplay between positive powers of and negative ones. In section 3.4 we have determined the reduced soft amplitude, thus we are in a position to explicitly write down the soft part of the infrared-renormalized amplitude at → 0, according to eq. (3.39a). Inserting eqs. (3.6) and (3.36) into eq. (3.39a) we get where we recall that x = L α s /π, and in the first line, corresponding to Z (−) NLL , we have dropped the −1 term in the numerator inside the square brackets, which does not generate any poles (see discussion following eq. (3.6)). This expression can be rearranged as follows: first of all, by collecting a factor e B 0 ( )−1 2 x C A we get We see at this point that the second line nicely cancel the poles from the first line. Furthermore, given that , it is safe to set to one all exponentials containing the factor 1 − B 0 ( ). We thus obtain NLL given in eq. (3.35). For later reference we also expand eq. (3.42) to the first few orders in perturbation theory, obtaining (recall It is interesting note that ζ n values with even n originate solely from the expansion of the factor R( ) in eq. (3.42), while the expansion of the factor∆ NLL generates only ζ n values with odd n. The latter property of∆ (+) NLL makes this function compatible with the class of zeta values we will encounter considering the two-dimensional amplitude in section 5.
In summary, according to (3.38) the infrared-renormalized amplitude is given as a sum of two terms: H s , computed in this section using the soft approximation, plus H h , which is identical to the hard part of the reduced amplitude (see eq. (3.39b)). The latter is infrared finite and originates in the hard wavefunction, which can be computed directly in two transverse dimensions. Let us turn now to evaluate it.

BFKL evolution in two transverse dimensions
As discussed in the introduction and in section 2, much of the complication of solving the BFKL evolution stems from the d-dimensionality of the Hamiltonian. Recalling that the two-reggeon wavefunction is finite at any loop order and that singularities are exclusively created by integration near the soft limit, it should be clear that no regularisation is required if we (a) only care about finite terms, and (b) remove any soft kinematics from the last integration. The latter condition is fulfilled by construction, having defined the split between the hard and soft wavefunctions (3.10) subject to the condition (3.11): the vanishing of Ω Our task in this section is therefore to compute Ω (2d) h (p, k). We do so by iteratively applying the Hamiltonian of eqs. (2.10) and (2.11), according to eq. (2.15). We keep the kinematics general, but in contrast to section 2, we work strictly in two transverse dimensions. To exploit the advantage of two-dimensional kinematics let us view the Euclidean momentum vectors k, k and p as complex numbers where the real and imaginary parts are the components of the corresponding momenta and introduce new variables z, w ∈ C according to Since the wavefunction is a function of Lorentz scalars (i.e. squares of momenta) it will be symmetric under the exchange z ↔z withz the complex conjugate of z. In particular, Ω h (p, k) depends on the two ratios These relations also clarify that the symmetry under interchanging the two Reggeons, eq. (2.24), corresponds to z → 1/z, and specifically, the two soft limits where one or the other Reggeon is soft correspond respectively to z → 0 and z → ∞. The limit z → 1 instead represents maximally hard kinematics, where both k 2 and (p − k) 2 are much larger than p 2 .
In the new variables the BFKL kernel (2.12) reads . (4.5) Furthermore, in the limit → 0, J(p, k) of eq. (2.13) becomes and the measure reads (4.7) Here, d 2 w ≡ dRe(w) dIm(w) where the real and imaginary part of w are to be integrated from −∞ to +∞, in accordance with eq. (4.2).
In applying BFKL evolution we employ the same notation as in the d-dimensional case but add the subscript "2d" to avoid confusion. In particular, from here on we express the two-dimensional hard wavefunction as Ω (2d) h (p, k) = Ω 2d (z,z). We expand it as in eq. (2.14), where we take B 0 (0) = 1, i.e.
where the coefficients of increasing orders are related by the action of the Hamiltonian according to eq. (2.15), which now reads: Plugging in the above expressions we find the two parts of the Hamiltonian to bê where In the next section we proceed to solve for the wavefunction Ω 2d by iterating the twodimensional Hamiltonian (4.9).

The two-dimensional wavefunction
It is useful to settle on a language before diving into the iteration of the two-dimensional wavefunction. To this end we introduce the class of iterated integrals dubbed single-valued harmonic polylogarithms (SVHPLs), which were first described by Brown in ref. [26]. Since then, several applications of SVHPLs in computing scattering amplitudes have been found, in particular in the context of the high-energy limit, e.g. [30][31][32][33][34][35], and in the context of infared singularities in general kinematics [36,37]. Here we will show that these functions also form a suitable basis for expressing the two-dimensional wavefunction Ω ( ) 2d (z,z) defined above.
As the name suggests, single-valued harmonic polylogarithms are single-valued functions which can be written as linear combinations of products of harmonic polylogarithms (HPLs) of z with HPLs ofz. We shall denote SVHPLs by L σ (z,z) where σ is a sequence of letters, typically zeros and ones. 3 The letters are said to form an alphabet, {0, 1}, and σ is, by analogy, referred to as a word. The length of a word is often called the (transcendental) weight of the SVHPL.
SVHPLs are the natural choice for the two-dimensional BFKL evolution, since j(z,z) of eq. (4.6) belongs to this class, 12) and the action of the HamiltonianĤ 2d,i preserves single-valuedness when acting on a singlevalued function. This can be expected on general grounds: any complex pair z,z identifies a point in the Euclidean transverse momentum plane. Physically there cannot be branch cuts in the Euclidean region; this, by definition, guarantees single-valued results. Indeed, singlevaluedness may be confirmed at every step of the iteration. Determining the wavefunction is greatly simplified by working directly with SVHPLs; we briefly summarise their main properties, which will be used below, in Appendix B.
As noted upon introducing the variables z andz in (4.2), the two-dimensional wavefunction is symmetric under z ↔z. In addition, as mentioned following (4.3), owing to the symmetry upon interchanging the two Reggeons in eq. (2.24), the wavefunction is invariant under simultaneously swapping z ↔ 1/z andz ↔ 1/z. Both these symmetries are easily verified by looking at eqs. (4.6) and (4.11a), where, for the latter symmetry, one changes the integration variables w → 1/w,w → 1/w. We will use these properties to simplify the iteration of the wavefunction as well as its results in section 4.3.
The evolution of the wavefunction in strictly two transverse dimensions according to (4.9) has the following basic characteristics. Firstly, iteratingĤ 2d,m amounts to multiplying by j(z,z) and therefore evaluating shuffle products of SVHPLs. Secondly, each application ofĤ 2d,i adds one layer of integration such that Ω can be written as a linear combination of SVHPLs of weight − 1. A method to calculate the convolution in eq. (4.11a) in terms of residues was described in chapter 6 of Ref. [35]. Here we develop an alternative method: we translate the action of the Hamiltonian into a set of differential equations, which we then solve in terms of SVHPLs.
Suppose we wish to compute the action of a linear operatorÔ, which may involve integration, on a function Ψ(z,z). Assume now that we find a differential operator ∆, which is linear in logarithmic derivatives with respect to z andz, with the following properties: ii. ∆Ψ is a pure function with a weight that is lower than Ψ by one unit. (4.13b) and we can computeÔ [Ψ(z,z)] by integrating the differential equation (4.14), assuming that the r.h.s. is known explicitly. If it is not the procedure can be applied recursively, i.e.
until the r.h.s. is simple enough to be calculated. After each integration a constant has to be fixed, e.g. by matching to known boundary conditions. Importantly, because ∆ is assumed to be linear in derivatives with respect to z andz, solving the differential equation amounts to computing a one-dimensional integral. This may be contrasted with the original integral in (4.11a) which is two-dimensional. Given (4.13b), solving this differential equation is straightforward, and the result remains in the class of HPLs, (see eq. (A.1)). The same applies for the class of SVHPLs: to solve the differential equation within this class, we first integrate its holomorphic part according to eq. (B.9), and subsequently recover the full result, depending on both z andz, by applying the single-valued map s defined in eq. (B.10). Having outlined the general approach let us see how it is implemented in practice to solve for the wavefunction in (4.11a).
Let us start by considering theÔ in eq. (4.14) to coincide with the two-dimensional HamiltonianĤ 2d,i (we will see below that the final procedure involves considering parts of the HamiltonianĤ 2d,i in turn). The most natural candidate for the operator ∆ in eq. (4.14) is ∆ 1 = z d dz , since condition (4.13a) is satisfied, as we now show. For generic values of w and z one finds using eq. (4.5) This implies that z d dz commutes with the Hamiltonian, fulfilling condition (4.13a). However, some extra caution is needed here: the complexconjugate pairs w,w and z,z cannot be treated as independent variables everywhere. Derivatives w.r.t. those variables receive additional contributions from the non-holomorphic or singular points of the function they act on. These "anomalies" are captured by the twodimensional Poisson equation namely, by contributions of the form d dw with c a complex number. The two-dimensional δ function in the above equations fixes both the real and the imaginary part of its argument such that for some function f , cf. the remark below eq. (4.7). For easy bookkeeping let us split a derivative into its regular part ("reg"), which is correct in the holomorphic regime, and its contact terms ("con"), governed by eq. (4.19). Eq. (4.16) therefore correctly reads which modifies eq. (4.17) to give We shall continue to refer to the behaviour in eq. (4.22) as the commutativity of z d dz and H 2d,i even though we implicitly mean commutativity modulo contact terms. Note, that the presence of the contact terms does not conflict with the strategy outlined above; each contact term contains a (two-dimensional) δ-function which makes the integral on the r.h.s. of eq. (4.22) easy to evaluate.
We will derive the explicit form of the contact terms towards the end of this section, at which point eq. (4.22) will become directly usable for determining the action ofĤ 2d,i on the wavefunction Ψ. Before doing that, however, we turn our attention to condition (4.13b). Concretely in eq. (4.22) the requirement is that z d dz Ψ should be a pure function of weight ones less than Ψ itself. We find that the operator z d dz , upon acting on any SVHPL of the form L 0,σ (z,z), does indeed yield such a pure function, so eq. (4.22) becomes: where we have used eq. (B.1). On the other hand, z d dz does not have the same effect when acting on an SVHPL L 1,σ (z,z), where one obtains instead which does not fulfil the condition (4.13b). One may be tempted to use (1 − z) d dz instead but, unfortunately, this operator does not commute withĤ 2d,i .
The solution is to first split the HamiltonianĤ 2d,i =Ĥ 2d,i 1 +Ĥ 2d,i 2 witĥ and where K 1 (w,w, z,z) + K 2 (w,w, z,z) = K(w,w, z,z), cf. eq. (4.5). This split is useful because it opens the possibility of identifying different differential operators ∆ i that commute with the separate components of the HamiltonianĤ 2d,i 1 andĤ 2d,i 2 , and yield a pure function when acting directly on L 0,σ (z,z) or on L 1,σ (z,z), thus simultaneously fulfilling both conditions in (4.13).
Regarding the commutation relations, condition (4.13a), it is straightforward to verify that the following four relations hold, up to contact terms: Let us therefore define the following three differential operators: and show that we can arrange the wavefunction, which is a linear combination of L 0,σ (z,z) and L 1,σ (z,z), such that condition (4.13b) would also be fulfilled.
To this end, let us first note that upon acting on L 0,σ (z,z) with either of the two parts of the Hamiltonian we have (using (4.27)): just as in (4.23). Thus, the remaining challenge is to handle terms containing L 1,σ (z,z); this is where the additional flexibility of splitting the Hamiltonian pays off. Let us consider first the simplest case ofĤ 2d,i 2 where we obtain NowĤ 2d,i 2 Ψ can be readily integrated for any Ψ using (4.29) and (4.30). Turning to considerĤ 2d,i 1 , let us write and use the linearity of the Hamiltonian to act with it on (L 1,σ + L 0,σ ) and (−L 0,σ ) separately. We may now apply respectively the differential operators ∆ 3 and ∆ 1 of (4.28) to these terms. With eq. (4.27a) and (B.1) one can easily verify that they produce the desired pure functions of lower weight in accordance with (4.13b): z(1 − z) d dz Ĥ 2d,i 1 (L 0,σ (z,z) + L 1,σ (z,z)) =Ĥ 2d,i 1 [L σ (z,z)] + (contact terms) (4.32) Using (4.32) along with (4.29) we see that alsoĤ 2d,i 1 Ψ can be integrated for any Ψ. Thus, by splitting the Hamiltonian and the wavefunction in a convenient way, we were able to identify linear differential operators that admit both requirements in (4.13).
In order to complete the process of setting up the differential equations let us now return to derive the explicit form of the contact terms. First, let us write eq. (4.22) for general ∆ i = f i (z) d dz and the two parts of the split Hamiltonian, where, according to eqs.  In computing the contact terms in (4.33) we note that the f i (z) (4.28) are functions of z only whilst being independent of the complex conjugatez. According to eq. (4.19) this implies that for n = 1, 2, and thus (4.33) becomes: Consequently, we only have to consider the following four derivatives, where in eqs. (4.37b) and (4.37d) we have dropped terms proportional to δ 2 (z), restricting our calculation to z = 0 (we emphasise that z is an external variable so this can be consistently done). Due to the sum of contact terms inside the curly brackets in eq. (4.36) the terms proportional to δ 2 (w − z) = δ 2 (z − w) in eqs. (4.37a)-(4.37d) cancel identically, so the remaining contact-term contributions are only at w = ∞ for K 1 and at w = 0 for K 2 .
Using the corresponding δ functions to turn the integrals over w in (4.36 ) into evaluation of limits at infinity and at zero respectively we finally obtain: (4.39) This equations will be used in the next section to determine the wavefunction.

Differential equations and an iterative solution for the wavefunction
Finding the differential equations is now simply a matter of compiling together the results of the previous section. Starting with the easiest case, ∆ 1Ĥ2d,i n L 0,σ , we notice that with f 1 (w) = w both the w → ∞ limit in eq. (4.38) and the w → 0 limit in eq. (4.39) vanish, and thus there are no contributions from contact terms in either of these cases. Dividing by f 1 (z) = z to arrive at Next consider the case ∆ 2Ĥ2d,i 2 L 1,σ , corresponding to eq. (4.30). Here f 2 (w) = 1 − w and eq. (4.39) yields where we have divided by f 2 (z) = 1 − z and used the shorthand [. . .] w,w→0 to denote the w,w → 0 limit of the functions inside the square brackets. This term can, in fact, be dropped as it always contains a single SVHPL whose indices feature (at least) one "1" and, thus, is equal to zero in the limit. The last case, ∆ iĤ2d,i 1 L 1,σ , is governed by eqs. (4.32) and (4.29), using the wavefunction split of eq. (4.31). Considering in turn the action of eq. (4.38) on (L 1,σ (z,z)+L 0,σ (z,z)) with f 3 (w) = w(1 − w) and on (−L 0,σ (z,z)) with f 1 (w) = w, we derive two separate equations, which we then combine using the linearity of operatorsĤ 2d,i 1 and d dz to obtain with [. . .] w,w→∞ the w,w → ∞ limit of the functions inside the square brackets. Taking this limit requires some careful analytic continuation of the relevant HPLs to ensure that w andw stay complex-conjugate as they approach infinity.
Because the HamiltonianĤ 2d,i and its componentsĤ 2d,i n are linear operators one can sum up the above equations (4.40)-(4.42) and recombineĤ 2d,i 1 +Ĥ 2d,i 2 →Ĥ 2d,i obtaining more compact expressions: These differential equations compactly represent the action of the HamiltonianĤ 2d,i according to eq. (4.11a). By solving them we are able to effectively bypass the computation of the two-dimensional integrals in the latter equation.
Since the differential equations only fix the z dependence of the (wave)function -which is a function of both z andz -a small detour is necessary to recover the action ofĤ 2d,i on SVHPLs: we take the holomorphic part of a given SVHPL, integrate it w.r.t. z according to the differential equations in (4.43), and then reconstruct the functional dependence onz by requiring the result be single-valued. This ultimately amounts to simply replacing For more details on this procedure see appendix B.1.
After each integration we need to fix an integration constant. We find that this is conveniently done by matching with the soft limit. Specifically, it is convenient to consider the soft limit where k 2 /p 2 = zz tend to zero. For small z,z, only SVHPLs with all-zero indices can give non-zero contributions; these correspond to powers of logarithms: In eq. (3.19) we calculated the action of the small-k (or soft) HamiltonianĤ s on powers of ξ = (k 2 /p 2 ) − . The action ofĤ i in the soft limit can be isolated by looking at the coefficient of 2C A − T 2 t and thus isĤ whereB m ( ) is given in eq. (3.20). Expanding both sides in and matching powers of δ = m in the limit → 0 lets us extract the action ofĤ i in the soft limit on any given power of log(k 2 /p 2 ) = log(zz). For reference, we find etc., from which we observe that the integration constants exhibit a very simple pattern. Specifically, they only contribute single (ordinary) zeta numbers because they are generated upon expandingB m ( ) which is a product of gamma functions.
We can now calculate the action ofĤ 2d,i on any SVHPL by iteratively solving the differential equations (4.43a) and (4.43b), starting from the lowest-weight functions, L 0 and L 1 . Effectively, we have set up an algorithm for calculating the two-dimensional wavefunction to any loop order. Due to the finiteness of the wavefunction it is straightforward to verify the results numerically: We integrate eq. (4.11a) numerically and compare to the analytical result for a number of randomly generated pairs z,z. Specifically, with w = w 1 + iw 2 and z = z 1 + iz 2 the action ofĤ 2d,i (4.11a) can be written where Ψ(z,z) is a (linear combination of) SVHPL(s). This type of integral is readily evaluated numerically in e.g. Mathematica.
Interestingly, a new type of transcendental number appears for the first time in the twelve-loop wavefunction -a so-called multiple zeta value (MZV). While it is no surprise that MZVs do not appear at lower loop orders as we explain in the following two paragraphs, the fact that they do appear starting at twelve loops is a non-trivial statement with numbertheoretical implications.
MZVs are the values of HPLs evaluated at special points, typically their branch points z = 1 or z → ∞, for example 4 H 0,0,0,0,1,0,0,1 (1) = H 5,3 (1) = ζ 5,3 . It turns out that SVHPLs only cover a subset of all MZVs when evaluated at z =z = 1 or z,z → ∞ and we refer to this subset as single-valued multiple zeta values. They are discussed in detail in refs. [27,29] where the authors show that, up to weight ten, the algebra of single-valued MZVs is generated by ordinary (odd) zeta numbers ζ n . At weight eleven, however, a new type of number appears, alongside the expected ζ 11 . We shall call it 5 g 5,3,3 and it is defined by where ζ 5,3,3 = H 5,3,3 (1).
There are two sources that contribute (multiple) zeta values to the wavefunction: the integration constants fixed by the soft limit and the w,w → ∞ limit in eq. (4.43b). The former are generated by expanding gamma functions, cf. eq. (4.46) with eq. (3.20), and can therefore contribute only single (ordinary) zeta numbers. The value of the large-w,w limit instead does generally involve (single-valued) multiple zeta values. We note that it is guaranteed to multiply the weight-one SVHPL L 1 (z,z) which is generated by the denominator, 1 − z, upon integrating the differential equation (4.43b). Being the sole source of (single-valued) MZVs, we conclude that such zeta values of weight w can only occur starting at the next loop order, i.e. = w + 1. Specifically, this explains why g 5,3,3 , which is weight 11, cannot appear at loop orders < 12. Indeed, we find that g 5,3,3 is accompanied by L 1 in the twelve-loop wavefunction: 4 MZVs use the collapsed notation, cf. eq. (A.4) in appendix A. 5 Brown [27] refers to it as ζsv(3, 5, 3) while Schnetz [29] calls it g335.
According to ref. [27] (cf. eq. (7.4) there) two more such numbers have to be introduced at weight 13 and, using the same logic, we anticipate that they make an appearance in the 14-loop wavefunction. Indeed, defining g 5,5,3 = 10ζ 2 2 ζ 9 + 275 2 ζ 2 ζ 11 + 5ζ 5 ζ 5,3 + ζ 5,5,3 The observed term g 5,3,3 L 1 (z,z) at twelve loops immediately rules out the possibility to find a closed-form expresson for the two-dimensional wavefunction in terms of gamma functions as was done in the soft limit (3.21). The non-zero coefficients of g 5,5,3 L 1 (z,z) and g 7,3,3 L 1 (z,z) at 14 loops may be seen as hint that indeed all single-valued MZVs appear in the two-dimensional wavefunction -when and as soon as the weight, i.e. loop order, allows for it. We will, in fact, encounter a contribution proportional to g 5,3,3 in the amplitude at eleven loops. We will thus return to discuss single-valued MZVs when interpreting our results for the amplitude in section 5.3.
Before we press ahead and compute the amplitude it is worthwhile exploring the aforementioned symmetries of the wavefunction in some more detail and we do so in the next subsection. This will ultimately lead to a better understanding of the iteration in two dimensions and enable us to calculate it to even higher loop orders.

Alphabets and symmetries
Throughout this paper we have tried to exploit the symmetries of the BFKL evolution to aide calculations and simplify expressions. In this section we explore to what extent symmetries can guide us in the two-dimensional limit. As mentioned in section 4.1, in two dimensions, the wavefunction is invariant under two transformations: complex conjugation and inversion of the arguments. The latter, i.e. the fact that Ω 2d (z,z) = Ω 2d (1/z, 1/z), corresponds to eq. (2.24), i.e. to swapping the two reggeons, and was used, for example, to identify the two soft limits in section 3. In the present context, it inspired us to introduce a new alphabet for SVHPLs, as we now explain. Instead of 0 and 1, corresponding to integration over d log z and d log(1 − z), respectively, we shall use a and s. They are associated with integration over d log z and d log z/(1 − z) 2 and thus behave antisymmetrically and symmetrically, respectively, under z → 1/z. In particular The leading-order wavefunction simplifies to Ω (1) 2d = 1 2 C 2 L s (z,z), and at higher orders, the z → 1/z symmetry implies that the antisymmetric letter a would only ever appear an even number of times.
Let us now consider the evolution directly in terms of this alphabet. Using the letters a and s simplifies j(z,z) = L s (z,z)/2 of eq. (4.6) and hence the action ofĤ 2d,m in eq. (4.11b), which now amounts to shuffling an s into the indices of the function it acts on (and multiplying by a 1 2 ), for examplê H 2d,m L a,s,a,s (z,z) = 1 2 L s,a,s,a,s (z,z) + L a,s,s,a,s (z,z) + L a,s,a,s,s (z,z). The action ofĤ 2d,i has a much richer and more complicated structure. However, we notice that at symbol level, i.e. keeping only the highest-weight SVHPLs, it simply amounts to replacing s → ss − aa and multiplying by − 1 4 , for example, H 2d,i L a,s,a,s (z,z) = − 1 4 L a,s,s,a,s (z,z) − L a,a,a,a,s (z,z) + L a,s,a,s,s (z,z) − L a,s,a,a,a (z,z) where Σ sub contains products of subleading-weight SVHPLs and zeta numbers, i.e. terms like L σ ζ n 1 · · · ζ nm with |σ| + n 1 + · · · + n m = 5 and |σ| < 5 in the above example. This replacement rule can be derived from the differential equations (4.43a) and (4.43b), as we now explain. To this end, let us consider the two casesĤ 2d,i L a,σ andĤ 2d,i L s,σ in turn. Considering the former, due to the equivalence of the letters 0 and a, eq. (4.43a) immediately gives the action on L a,σ d dzĤ 2d,i L a,σ (z,z) =Ĥ 2d,i L σ (z,z) z . (4.59) The simple recursive nature of this equation implies thatĤ 2d,i does not affect the a indices of a SVHPL and can, at most, generate subleading terms Σ sub through integration constants, cf. eqs. (4.47a)-(4.47e).
The action on L s,σ can be written as where at each step we have used Σ sub to collect subleading terms into. The first term in the final expression is again an inert term, like the one encountered in eq. (4.59). The following term however, creates two leading-weight terms which, upon integration, yield − 1 4 (L s,s,σ − L a,a,σ ) and hence confirm the pattern described above eq. (4.58). Note that by the recursive nature of the differential equation this applies (separately) to every letter s in the word (s, σ), not just the first one (see e.g. eq. (4.58)).
In the following we show that it is possible to unravel the recursive definition ofĤ 2d,i beyond symbol level. The Σ sub terms in the above equations are generated by two independent and additive sources: the w,w → ∞ limit in eq. (4.43b) and the constants of integration as shown in eqs. (4.47a)-(4.47e). Let us denote them Σ sub(∞) and Σ sub(0) , respectively, with their sum equalling Σ sub . Empirically we observe that Σ sub(0) follows a simple pattern when using the {a, s} alphabet: (4.61) with Σ lead now the leading-weight SVHPLs governed by eq. (4.60). Σ sub(∞) in turn can be summarised bŷ × L a,w j+1 ,...,w −1 (z,z) + L s,w j+1 ,...,w −1 (z,z) z,z→∞ . (4.62) In both these equations the final term in the sum needs to be interpreted with care: in eq. (4.61), for j = one obtains L w 1 ,...,w 0 ≡ 1 and in eq. (4.62) for j = − 1 one obtains in the second factor L a,w ,...,w −1 (z,z) + L s,w ,...,w −1 (z,z) ≡ L a (z,z) + L s (z,z). Observe that w j ! = s in eq. (4.62) is a necessary yet not sufficient requirement for a non-zero contribution. Being based on observations, the patterns described in eqs. (4.61) and (4.62) need to be verified against the wavefunctions computed in the previous section. We find perfect agreement with the wavefunction up to and including 13 loops, and are thus confident that the above description is correct.
By introducing the {a, s} alphabet we have accounted for the symmetry of the wavefunction under inversion, z → 1/z, at symbol level, i.e. as far as leading-weight terms are concerned. Our basis of SVHPLs respects neither this nor the invariance under complex conjugation at function level: in general L σ (z,z) = L σ (1/z, 1/z) and L σ (z,z) = L σ (z, z). Expecting further simplifications we will therefore construct a set of symmetrised functions in the remainder of this section.
In the following we heavily use relations between SVHPLs under a standard set of variable transformations. We summarise the most important aspects of these relations in appendix B.2. Quintessentially, these relations determine the coefficients c w in L σ (g(z), g(z)) = w c w L w (z,z) where the sum runs over all words up to weight |σ| and, in the present case, g(x) = 1/x or g(x) =x.
Let us define with σ a word belonging to an alphabet of one's choosing. In the following we stick with the {a, s} alphabet. We stress that the set of Fs does not span the space of SVHPLs but it does cover the entire space of wavefunctions. Due to the symmetries of the wavefunction and thus Ω 2d (z,z) = 1 4 (Ω 2d (z,z) + Ω 2d (z, z) + Ω 2d (1/z, 1/z) + Ω 2d (1/z, 1/z)) (4.65) one can simply replace L σ (z,z) → F σ (z,z) to go from the L to the F basis. It may therefore not be immediately obvious how eq. (4.63) simplifies the results. Indeed, it requires a few more steps to showcase the advantages of a symmetrised basis. Firstly, the wavefunction in the L basis contains functions whose indices feature an odd number of the letter a. Their leading-weight components are antisymmetric under z → 1/z because d log z = −d log 1/z (4.66) Converted to F functions they are hence zero at symbol level or, in other words, equal to products of lower-weight SVHPLs and zeta numbers. This can be turned into a recursive algorithm that successively removes all odd-a functions. Schematically, 1. Consider the wavefunction at a given order and replace L σ (z,z) → F σ (z,z) 2. Choose an F σ (z,z) where σ contains an odd number of a letters. Plug in definition (4.63) and rewrite SVHPLs as functions of z,z using the rules in appendix B.2. The resulting SVHPLs will be of lower weight than the original F σ , multiplied by zeta numbers.

Repeat steps 2 & 3 until a fixed point is reached and only functions with an even number of a letters remain.
Note that step 3 is valid for the same reason it was legitimate to replace L σ (z,z) → F σ (z,z) in the wavefunction, cf. eqs.  Secondly, we may combine F σ (z,z) and Fσ(z,z) withσ the word σ reversed, at the cost of generating subleading terms. This is due to the following identity of SVHPLs: For a function F σ this entails due to the invariance under complex conjugation. Besides removing nearly half of the F functions we find the generated subleading terms to sometimes reduce but never increase the complexity of a given expression. For the procedure to be algorithmic one chooses which letter to cumulate in the left (or right) half of a word. For the wavefunction up to four loops and with the same abbreviations as in eqs. (4.49a)-(4.49d) we find where we used eq. (4.69) in favour of words that start rather then end with the letter s. Further results up to weight 13 can be found in the ancillary file 2Reggeon-wavefunction-Fsa-Basis.txt. Indeed, comparing the results in eqs. (4.70a)-(4.70d) to the wavefunction in terms of standard SVHPLs (and the standard {0, 1} alphabet) in eqs. (4.49a)-(4.49d) shows the benefits of the new basis. In terms of F functions the wavefunction takes not only a very compact form and is expressed in terms of fewer functions, it also removes subleading terms in some cases, like the − 3 16 L 1 ζ 3 in the coefficient of C 3 1 C 2 at four loops (4.49d).

Finite corrections to the amplitude from two-dimensional evolution
We now have an algorithm for the calculation of the wavefunction Ω 2d to any loop order, and we shall use it for the computation of the amplitude. Let us recall from section 3 that the soft part has been fully determined, and our goal here is the calculation of the hard part of the amplitude, as defined in eq. (3.14b). This, in turn, requires the hard part of the two-dimensional wavefunction, which according to eq. (3.12) is obtained by subtracting the d = 2 limit of the soft wavefunction from the full (two-dimensional) wavefunction Ω 2d of the previous section. To this end we first define where taking the limit simply corresponds to selecting the leading O( 0 ) terms in Ω s (p, k). Notice that within the d = 2 limit we switch to the two-dimensional variables z andz of eq. (4.2), and the single-valued logarithm L s (z,z) = log zz (1−z) 2 (1−z) 2 defined in eq. (4.56). Having used the symmetrised soft wavefunction we land directly in the class of SVHPLs used in the two-dimensional computation of section 4.2, and in this way the computations of Ω 2d and Ω (2d) s are entirely compatible. We note that with the replacement in (5.1) the two-dimensional soft wavefunction in (3.23) becomes a polynomial in L s (z,z) at any given order. According to eq. (3.24), these terms exponentiate and can be resummed to all-orders. Upon applying the changed of variable of (5.1) this resummed expression is Using the results in section 4 for Ω 2d and the expansion of eq. (5.2) for Ω , and we can proceed to determine the hard part of the amplitude order by order, according to eq. (3.14b). To this end, recall that the hard wavefunction Ω (2d) h is guaranteed to integrate to finite terms only, hence it can be integrated in strictly two dimensions. Applying the limit → 0 in eq. (3.14b) to the integrand and the integration measure using the variables z andz (cf. eqs. (4.3) and (4.7)) we obtain: where, in practice, we loop-expand both the wavefunction and amplitude, as was done in eq. (2.18). The next two subsections are dedicated to the computation of the integral in eq. (5.3), thus determining the hard component of the reduced amplitude order by order. In section 5.3 we combine the soft and hard components of the reduced amplitude according to eq. (3.13), and finally in section 5.4 we similarly combine the soft and hard components of the infrared-renormalized amplitude using eqs. (3.38) and (3.39).
To set up the computation of eq. (5.3) let us define and introduce in turn two independent methods for computing these integrals. For the sake of simplicity of notation, given that the entire computation is done in two dimensions, we shall now drop the (2d) superscript, and refer to the integrand in (5.4) as Ω h (z,z). Similarly, while (5.4) is applied order-by-order, in describing the methods we refrain from using an index for the loop order on either side of (5.4). The first method, described in section 5.1 below, is based on using the known analytic structure of the wavefunction, in order to convert the two-dimensional integral into an integral over the discontinuity of the wavefunction. It was inspired by the calculations described in section 7.1 of ref. [50]. The second method, presented in section 5.2 below, relies on the symmetry of the wavefunction under inversion, z → 1/z,z → 1/z, and the action ofĤ 2d,i at fixed external points.

Method I: final integration using the discontinuity of the wavefunction
Let us define a regularised version of the integral I in eq. (5.4): where the cutoff δ is assumed to be small, δ 1. The introduction of δ may seem superfluous at this point as lim z,z→0 Ω 2d = lim z,z→∞ Ω 2d = lim →0 Ω s and thus, using eq. (5.2), lim z,z→0 Ω h = lim z,z→∞ Ω h = 0; more precisely, Ω h vanishes linearly in zz in the soft limit, up to logarithms, rendering the integral in (5.4) convergent, and the difference I − I reg = O(δ 2 ) (up to logarithms). The necessity of this cutoff despite this good convergence will become clear shortly.
The exclusion of the points {0, ∞} in (5.5) enables us to introduce polar coordinates such that zz = r 2 and z z = e 2iθ , as now all points in the integration region have a nonvanishing Jacobian: To proceed we express the angular integral in the latter as an integration in the complex y plane where y ≡ e iθ , getting y Ω h (ry, r/y) , where the contour runs along the unit circle. The method outlined in the following is based on deforming the contour in the complex y plane. Essential to this is the fact that the integrand, at any order, is expressed in terms of SVHPLs, whose analytic structure is well understood. These functions are single-valued as long as their arguments are complex conjugates of one another, namely as long as the contour in eq. (5.7) runs along the unit circle. Outside of this region, i.e. upon deformation the contour, the HPLs in Ω h (z,z) exhibit branch cuts where z ∈ [1, ∞] andz ∈ [1, ∞]. In the r, y coordinates of eq. (5.7) they correspond to cuts along the real axis in the complex y plane where y ∈ [1/r, ∞] and y ∈ [0, r], respectively.

ℜy
ℑy ℜy ℑy r 1 r Figure 6. Position of the branch cuts in z andz in the complex y-plane for r < 1 (left). The contour along the unit circle in eq. (5.7) can be deformed and, consequently, identified with the integral of thez-discontinuity (right), as written in eq. (5.8).
For r < 1 there is a branch cut-free interval (r, 1/r) through which the contour along the unit circle passes, cf. the l.h.s. of figure 6. The contour can consequently be shrunk until it corresponds to integrating thez-discontinuity of the wavefunction over y from 0 to r, cf. the r.h.s. of figure 6. We can now understand why it is necessary to work with the regularised integral I reg of eq. (5.5) instead of the original I of eq. (5.4): while the hard wavefunction Ω h (z,z) vanishes at 0 and ∞, its discontinuity, in general, does not. In other words, the contour deformation introduces spurious divergent terms and the cutoff introduced in eq. (5.5) regularises them. .
For r > 1 the branch cuts of z andz overlap. However, the discontinuity cancels identically in the interval (1/r, r). Repeating the procedure, we again identify the contour integration with integrating thez-discontinuity of Ω h (z,z) over y, this time, from 0 to 1/r.
In total, having modified the contour in (5.7) we find where the two terms correspond respectively to r < 1 and r > 1. It is clear from the outset that they are equal: this corresponds to splitting (5.5) at zz = 1, which admits (z,z) ↔ (1/z, 1/z) symmetry. In the second line of (5.8) we reverted to the variablē z = r/y in the first integral and z = ry in the second; in the third we changed the order of integration before reverting to z = r 2 /z in the first integral andz = r 2 /z in the second; finally in the last line we defined x = 1/z in the both integrals. Let us now discuss the evaluation of the final expression in eq. (5.8), where the integration region of two terms is depicted as the white area in figure 7. In order to perform the integration it is useful to view the integrals (cf. the r.h.s. of figure 7) as the integral over a square plus (the integrals over) two wedges (5.10) minus two small triangles (5.11) where both z and x are small. Next we would like to evaluate each of these contributions, distinguishing between finite, δ-independent terms, and logarithmically divergent cut-off dependent terms.
The discontinuity w.r.t.z of Ω h (z, 1/x) evaluates to HPLs of z and x. I A (δ) of eq. (5.9) thus immediately evaluates to HPLs at 1, giving rise to MZVs, and at δ 2 ; the latter contain logarithmically divergent terms in δ. The first (second) integral in the expression of I B (δ) in (5.10) is calculated close to z = 0 (x = 0), cf. figure 7. One can therefore expand the discontinuity function in the integrand and discard terms suppressed by powers of z (x) keeping only powers of log z (log x). The inner integrals then yield powers of log δ 2 , log δ 2 x = log x + log δ 2 and log δ 2 z = log z + log δ 2 , respectively. The outer integrals thereupon generate MZVs from their upper limits; in addition it contains logarithmically divergent terms in δ. Contributions from the lower integration limits are dropped according to the (standard) regularisation of HPLs: A similar analysis of I C (δ) in eq. (5.11) reveals that only powers of log δ 2 are generated by the integrations over the two small triangles in figure 7.
Since the original integral I of eq. (5.4) is finite and I reg → I for δ → 0 all terms proportional to log δ 2 have to cancel between the three contributions I A (δ), I B (δ) and I C (δ). This enables us to derive a simplified integral in which the logarithmically divergent terms are absent altogether whilst giving the same finite terms: where all integrals are regulated according to eq. (5.12) and discz =1 [Ω h (z, 1/x)] z 1 and discz =1 [Ω h (z, 1/x)] x 1 refer to the aforementioned expansion of the integrand around small z and x, respectively. The first integral in eq. (5.13) reproduces all finite, cut-off independent terms in I A (δ) of (5.9), while the second and third ones reproduce, respectively, the finite terms in the two integral in I B (δ) in (5.10); finally, given that no cut-off independent terms are produce by I C (δ), it has no trace in (5.13). The above calculation is biased towards the discontinuity with respect toz which is purely a matter of choice. A similar calculation can be performed to get an answer in terms of the discontinuity with respect to z or a mixed expression that features both discontinuities.
This integration method was further checked as follows. Given a wavefunction (or SVHPL) we expand around z =z = 0 and change variables to the polar coordinates introduced above in eq. (5.7). The result is a sum of terms of the form r a y b log c (r 2 ) with rational constant coefficients and a, c ≥ 0 and b are integer powers. Integrating the azimuth over [0, 2π] then removes all terms that explicitely depend on y, i.e. that have b = 0. Next, we determine the rational coefficients in terms of harmonic numbers 6 . This enables us to perform the sum ad infinitum after we integrate term-by-term with respect to r.

Method II: final integration as an action of the Hamiltonian
The previous method, albeit straightforward on paper, is computationally demanding at high loop orders as it requires extensive use of analytic continuations of HPLs to calculate discontinuities. It turns out there is an easier way to perform the final integration, which lets us make use of our knowledge about the action of the Hamiltonian, established upon computing the wavefunction in section 4.
Consider the action ofĤ 2d,i (4.11a) on the wavefunction and set z =z = 1 under the integral.
Using Ω h (0, 0) = 0 one gets on the right-hand side: cf. eq. (4.5). It thus follows that (5.14), taken in the limit z,z → 1, yields: where in the second line we changed the integration variables in the first two terms -in the first using w → w/(w − 1), and in the second using w → 1 − w, and then factored out a common denominator. Given that the wavefunction is symmetric under inversion, Ω h (1/w, 1/w) = Ω h (w,w), the first and third terms in the last equation cancel and we find which can be readily identified with the integral in eq. (5.4) which we are interested to compute. We thus conclude that the integral in eq. (5.3), representing the hard wavefunction contribution to the reduced amplitude, integrated in exactly two dimensions, may be calculated with the methods we developed for the computation of the two-dimensional wavefunction itself, described in section 4.1. In practice one rewrites the wavefunction Ω h (1 − z, 1 −z) in terms of SVHPLs of z andz, then applies the Hamiltonian by solving the corresponding differential equations, and finally evaluates the resulting expression at z,z = 1. The last step produces the anticipated MZVs.
Method I, described in section 5.1, and method II outlined here show perfect agreement when applied to the wavefunction. However, we emphasise that while former may be applied on individual SVHPLs, the latter can only be applied to expressions which are symmetric under inversion of their arguments, cf. eqs. (5.17) and (5.18).

Results for the reduced amplitude
With the methods described in the previous sections it is straightforward to integrate the two-dimensional wavefunction and thereby compute the hard contribution to the amplitude, namely the finite terms not captured by the soft limit.
Before presenting our results let us recall the number-theoretic observations we made about the amplitude at the end of section 2. There, we claimed that the -loop amplitude (divided by B 0 , (2.4)) has two important number-theoretic properties: all of its terms have weight and there are no terms proportional to ζ 2 . We proved this statement for contributions from the soft limit in section 3, see below eq. (3.27h). We now show that it holds also for the hard contributions.
We begin by noting that the integrand in (5.4) is expressed as a pure function of uniform weight, written as sums of products of HPLs. We note that both methods for the last integral, in sections 5.1 and 5.2, increase the weight of the functions they act on by one, before evaluating the result at z =z = 1. In method I the action of the discontinuity first lowers the weight of its argument by one; this is then compensated by two consecutive integrations of a d log form, each raising the weight by one. Method II, in turn, applies the HamiltonianĤ 2d,i on the wavefunction after a variable transformation z → 1 − z. Changing the variables of an SVHPL does obviously not change its weight and the action of the Hamiltonian corresponds to integrating a first-order differential equation, which raises the weight of the operand by one.
SVHPLs at z =z = 1 evaluate to multiple zeta values (MZVs) of the same weight, as discussed following eqs. (4.49a)-(4.49d). We remind the reader that the ( − 1)-loop wavefunction consists of weight-( −1) SVHPLs and weight-( −1) products of SVHPLs and zeta numbers and conclude that the hard contributions to the -loop amplitude therefore have uniform weight . The absence of ζ 2 in the hard componentM (+, ) NLL,h is readily explained by the fact that SVHPLs can, by construction, only ever evaluate to odd zeta numbers, for any argument.
We start the discussion of the results by presenting the contributions that originate in the hard region. They are the immediate result of the previous sections and, through eight loops, read where we again used the shorthand notation for the colour factors, C 1 = 2C A − T 2 t and C 2 = C A − T 2 t . One may observe the aforementioned homogeneous weight property and absence of even zeta numbers. In fact, considering the first eight loop orders, one may get the false impression that each order contains just a single (product) of zeta numbers and that they are all single (ordinary) zeta numbers. Both these features are artefacts of looking at low weights, and a much richer structure will be revealed at higher loop orders, as we discuss shortly.
Given the identity in eq.  .43a). This will be done in the section 5.4 below. Before doing this let us combine the hard and soft components for the reduced amplitude itself, and comment further on some number-theoretic properties, as promised.
According to eqs. (3.13) and (3.14), the expressions for the full reduced amplitude through O( 0 ) can be easily obtained order-by-order summing the results for the soft amplitude provided in eqs. (3.27a)-(3.27h) and those for the hard amplitude in eqs. (5.19a)-(5.19h) above, where the former accounts for all infrared singularities plus some finite terms, and the latter for the remaining finite contributions. We obtain These result reproduce the one-to four-loop results of ref. [23] and of our numericallydetermined five-loop result in eq. (2.33). In the ancillary file NLL-reduced-amplitude.txt we provide the result for the soft, the hard and the full reduced amplitude up to 13 loops. Furthermore, the amplitude can now be calculated to any number of loops with the methods presented in sections 3, 4 and 5.
Similarly to the wavefunction at twelve loops (and above), the hard contributions to the amplitude (and thus the full amplitude itself) cannot be expressed in terms of ordinary zeta numbers beyond a certain loop order. In fact, most of what we discussed in the context of the wavefunction below eqs. (4.49a)-(4.49d) applies to the amplitude as well: Either of the two methods presented in section 5.1 and 5.2 requires us to evaluate SVHPLs at z =z = 1 and we hence expect the presence of (single-valued) MZVs starting from weight eleven. Indeed, the eleven-loop amplitude features a term proportional to g 5,3,3 , defined in eq. (4.50): Of course this term is entirely due to the hard component of the amplitude, as the soft one consists exclusively of ordinary zeta values (non-single-valued ones), as discussed in section 3. At twelve loops the reduced amplitude is again comprised of ordinary zeta numbers ζ n , as there are no weight 12 single-valued MZVs. Such numbers appear then again in the thirteen loop amplitude: where the single-values zeta numbers g 5,5,3 and g 7,3,3 have been defined in eqs. (4.52) and (4.53).
The fact that the MZV terms in eqs. (5.21), (5.22) and (5.23) appear in the elevenand thirteen-loop amplitude already excludes that there could be a simple all-order formula in terms of gamma functions for the reduced amplitude. This stands in sharp contrast to the contributions associated with the soft limit, both singular and finite, which can be resummed to all orders by means of gamma functions, as we have seen in section 3.

The infrared-renormalized amplitude
We conclude this section by discussing the perhaps most physically relevant infraredrenormalized amplitude (or hard function), which according to eq. (3.38), is obtained by summing the soft component, given to all orders by the closed expression in eq. (3.42), and the hard component, which according to eq. (3.39b), coincides with the hard component of the reduced amplitude,M (+) NLL,h . The latter can be determined to any loop order by following the methods discussed in sections 5.1 and 5.2, however a closed-form expression cannot be obtained as in case of the soft part of the infrared-renormalized amplitude. Thus, in practice we limit ourselves to determine this amplitude to 13 loops, and the result is provided in the ancillary file NLL-IR-renormalised-amplitude.txt. Here we provide a sample of the result (with H defined in eq. (3.1) and loop-expanded following eq. (2.17)), up to eight loops: (5.24h) NLL,h , the same number-theory properties discussed at the end of section 5.3 apply to the infrared-renormalized amplitude as well. In particular, resummation in terms of gamma functions is excluded.

Numerical analysis and convergence properties
The calculation developed in sections 3, 4 and 5 has allowed us to determine symmetries and general features of the wavefunction and the amplitude, as well as their exact analytic structure to fourteen and thirteen orders in perturbation theory respectively. We are now interested to perform a numerical analysis, and focus on features which are not directly evident from the analytic expressions, such as the qualitative behaviour of the wavefunction, the relative size of the soft and hard contributions to the wavefunction, and the convergence properties of the infrared-renormalized amplitude as an expansion in x ≡ L α s /π.

Wavefunction
Let us begin by analysing the wavefunction. Given its finiteness, we consider here the leading term in the expansion, i.e. the two-dimensional soft, hard and full wavefunctions, defined respectively in eqs. (5.1), (3.12) and (4.9).
As an example, in figures 8 and 9 we plot the coefficients Ω , Ω ( ) 2d , at third, fourth and fifth order in perturbation theory. In these plots we fix N c = 3 and consider specific colour representations, namely the singlet and the 27 representation, such that the Casimir operator in the t-channel evaluates to singlet :  We plot the wavefunction in the complex z plane, forz = z * . We observe that the soft and full wavefunctions exhibit peaks at z = 0; these are associated with the soft limit. Of course, by the z ↔ 1/z symmetry discussed in sections 4.1 and 4.3, there is an identical Figure 8. Soft, hard and full wavefunction in the complex plane Re(z), Im(z). Here we plot the singlet component of the wavefunction. The soft and the full wavefunction exhibit singularities at z = 0 and z = ∞, due to the z ↔ 1/z symmetry, (the latter is not visible in the plots). In addition, there is a singularity at z = 1, which appears also in the hard part of the wavefunction. Notice that the singularities at z = 1 partly cancel between the soft and hard wavefunctions, such that the full wavefunction exhibits a peak near z = 1 which is markedly smaller relative to the two separate contributions.
singularity at z = ∞ (which is not visible in the patch of the complex plane shown in the plot). In the way we separated between the soft and hard components, the two-dimensional hard wavefunction is strictly zero in these soft limits (see the discussion following eq. (5.5)).
All components of the wavefunction have singularities at z = 1. The z = 1 singularity represents rather different physics, where both Reggeons are hard, namely k 2 , (p−k) 2 p 2 . It is interesting to note that the singularity at z = 1 is always of opposite sign between the soft and hard wavefunction, such that these contributions cancel to a large extent in the full wavefunction. This observation allows us to conclude already that the soft approximation, although convenient for calculation purposes, does not provide a good numerical approximation for the full wavefunction away from the soft limit.
Focusing now on the full wavefunction, the singular behaviour near z = 0 and z = 1 at loop order can be described respectively by the leading logarithms in the two limits, ∼ c log (zz) and ∼ c log 1/(1 − z) 2 (1 −z) 2 , where both the magnitude and the sign of the coefficients c and c depends on the colour representation considered. Concerning the limit z = 0, the asymptotic behaviour is entirely determined by the soft wavefunction, given that Ω 2d h (0, 0) = 0. We obtain the coefficients c expanding the soft function in eq. (5.2) (compare eqs. (3.24) and (3.22)). Taking into account eq. (2.14), we find Given that (C A − T 2 t ) = C A = 3 for the singlet, while (C A − T 2 t ) = −C A − 2 = −5 for the 27 representation, this explains the sign-oscillating behaviour of the wavefunction near z =z = 0 for the singlet, and the constant sign of the 27 representation, which can be seen also in figures 8 and 9.
Determining the coefficients c is less trivial, given that near z = 1 also Ω 2d h contributes. An analysis of the asymptotic behaviour up to the 14th order allows us to deduce the pattern and extrapolate. We find: Once again, we see that the series has alternating or constant signs depending on the color representation. Specifically, it is sign alternating for the 27 representation, and it has constant sign for the singlet. Notice that both asymptotic expansion of the wavefunction near z = 0 and z = 1 can be summed using eq. (2.14). We obtain if T 2 t = 0.
where 0 F 1 and 1 F 1 are the confluent and Kummer's confluent hypergeometric function. These resummed expressions are valid only in the leading logarithmic approximation in zz and (1−z)(1−z), respectively. The generalization of (6.4) to include subleading logarithms of zz has been given in (3.24), while a closed form generalization of (6.5) is yet unknown.

Convergence of the loop expansion of the infrared-renormalized amplitude
Having computed finite contributions to the imaginary part of the amplitude to high loop orders we are in a position to investigate a very interesting theoretical question, namely the convergence properties of the perturbative expansion. Of course, this is done here at a fixedlogarithmic accuracy, namely considering the amplitude as a function of x ≡ L α s /π. The high-energy limit adds an interesting twist to the question of convergence, since within the a priori "perturbative regime" where α s (µ 2 ) is small (recall that µ 2 is naturally determined by the momentum transfer −t, and we assume s −t Λ 2 QCD ) high-energy logarithms L ∼ log |s/t| can be arbitrarily large, but then the effective expansion parameter x ≡ L α s /π becomes large. Thus, while there is no obvious reason why perturbation theory should break down, the question arises whether we can extend the validity of the calculation to large values of the expansion parameter x. In [25] we studied the infrared-divergent part of the amplitude in detail, and proved that these corrections exponentiate in terms of the soft anomalous dimension. We determined the latter to all orders in perturbation theory, and shown that it is an entire function, having an infinite radius of convergence in x.
We are now in a position to study the convergence of the infrared-renormalized amplitude H (+) NLL , which we determined analytically to the 13th order in section 5. For convenience, we introduce the amplitude Ξ and its coefficients Ξ ( ) , defined through  such that H (+, ) cf. eqs. (2.17) and (2.18). In the following we will use equivalent definitions also for the soft and hard parts of the infrared-renormalized amplitude coefficients. Numerical expressions for the coefficients of the infrared-renormalized amplitude Ξ (+, ) NLL up to thirteen loops can be obtained starting from the analytic expressions given in the ancillary files 7 , using the relations in eq. (6.1), and converting the multiple zeta values there into decimal numbers. We arrive at  for its hard part. We plot the partial sums of the soft, hard and full infrared-renormalized amplitude as a function of x respectively in figures 10, 11, and 12.
Considering the all-order resummed expression for the soft component of the infraredrenormalized amplitude in eq. (3.42), we can immediately conclude that it exhibits a finite radius of convergence. The radius of convergence can be identified as the position of the pole closest to the origin in the complex x plane, which we denote in what follows R. Inspecting eq. (3.42), and in particular the explicit expression for∆ (+) NLL in eq. (3.35), we see that the soft part of the infrared-renormalized amplitude has poles when the argument of the gamma functions in the numerator equals zero or negative integers. In general poles appear for both positive and negative x: this is an important point we shall return to below. The pole closest to the origin is determined by 1 − (C A − T 2 t ) x = 0, which in turn determines the radius of the convergence of the soft part of the infrared-renormalized amplitude to be R s = 1/(C A − T 2 t ) (the subscript "s" refers to the soft part of the infraredrenormalized amplitude). This corresponds to R s = 1/3 0.333 for the colour singlet infrared-renormalized amplitude, and R s = −1/5 = −0.2 for one in the 27 representation. The qualitative picture of convergence of the partial sums as a function of for any x < R s , and divergence beyond that point, can indeed be confirmed upon inspecting figure 10.
For the hard contribution to the infrared-renormalized amplitude, and thus also for the complete one, we do not have an all-order expression. Nevertheless, information on the radius of convergence can be extracted from the perturbative expansion by constructing Padé approximants of the infrared-renormalized amplitude. More specifically, we may use the partial sum of the infrared-renormalized amplitude at any order to construct a rational function of x, which reproduces the partial sum upon expansion. Here we choose to use Padé approximants of the form 8 : (6.14) With this definition the Padé approximant has two poles at x ± = −b 1 ± b 2 1 − 4b 2 /(2b 2 ). The pole closest to the origin provides a prediction for the radius of convergence of the series: R = min{x − , x + }. Of course, this prediction is expected to be reliable only upon considering sufficiently high orders, where the series approaches its asymptotic regime. The stability of the deduced value for the radius of convergence with respect to the order provides an indication of whether the asymptotic regime is reached.
Before describing the results a further comment is due regarding the sign of R. Strictly speaking, the radius of convergence would be the absolute value of R. Here, however, we are interested in keeping track also of the sign of the nearest pole, indicating whether the series is (asymptotically) of constant signs, or oscillating. We shall see that both scenarios are realised.
We test this method on the soft part of the infrared-renormalized amplitude, for which the all-order result is known, as discussed above. The results for the singlet and the 27 representation for N c = 3 are shown in Table 1. We see that in both cases the pole closest to x = 0 (x [1] s,− = 0.333 for the singlet and x [27] s,+ = −0.200 for the 27 representation) approximates very well the exact radius of convergence, R s = 1/(C A − T 2 t ).
singlet there is a highly stable nearest pole at x [1] h,− = 0.333. For the 27 representation, in turn, the stable pole at x [27] h,+ −0.19 is not always the one closest to the origin, due to the wide fluctuations of x [27] h,− . Finally, for the complete infrared-renormalized amplitude we summarise the results in table 3. Here we find highly stable results: x We conclude that Padé approximants based on partial sums of order = 10 through 13, yield fairly stable predictions for the poles. Naturally, ones still finds some fluctuations, which can be attributed to subasymptotic effects, but an overall consistent picture emerges, and we can deduce an approximate radius of convergence in each case from the position of the poles.  Table 3. Table summarising the values of x = αs π L at the poles of the Padé Approximants in eq. (6.14), considering the full amplitude at orders = 10 through 13 for the singlet and the 27 representation.
The final results of this analysis are summarised in table 4, where we compare the results for the soft part of the infrared-renormalized amplitude, deduced from the resummed result (which are highly consistent with the Padé approach), with those for the hard component and complete infrared-renormalized amplitude, which are both based solely on the Padé analysis. In the table we also provide an interpretation of the radius of convergence for the full infrared-renormalized amplitude in terms of the analytic dependence on the colour factors C 1 and C 2 ; this will be explained below.  Table 4. Summary table for the radius of convergence R of the expansion of the infraredrenormalized amplitude in powers of x = αs π L, determined by identifying the pole closest to x = 0 using Padé approximants. We use the shorthand notation The numerical results in the table indicate that the radius of convergence of the full infrared-renormalized amplitude is larger compared to both its soft and hard components. Indeed, better convergence is clearly observed looking at successive orders in the full hard function in figure 12 compared to its soft and hard components in figures 10 and 11, respectively. The interpretation is clear: the pole that limits the convergence of the soft component of the infrared-renormalized amplitude in the resummed expression, eq. (3.42), exactly cancels against a similar divergence in the hard component, hence the similar values of R for the soft and hard components in table 4. Upon cancelling the leading divergence, a subleading pole is exposed, which becomes the dominant obstruction for convergence of the full amplitude. This is of course another indication that the separation of the finite O( 0 ) terms between the soft and hard regimes is arbitrary; we have already seen that the soft wavefunction cannot approximate the full one away from the soft limit in figures 8 and 9.
Even more interesting is the observation that the sign of the first pole, R, which indicates whether the series is asymptotically sign-oscillating (R < 0) or of constant signs (R > 0), is negative for the full infrared-renormalized amplitude, while it may be either positive or negative for the separate soft and hard components, as can be seen in table 4. Upon resumming the perturbative expansion of the full infrared-renormalized amplitude, one expects a smooth extrapolation to high energies when taking the centre-of-mass energy large compared to the momentum transfer, s −t Λ 2 QCD . Given that the expansion parameter, x = αs(−t) π log s −t , gets large (and positive) in this limit, smooth extrapolation (of the resummed expression) to high energies can only be consistent with a finite radius of convergence if the series is sign-oscillating, or put in stronger terms: if all the singularities of the resummed infrared-renormalized amplitude are locate away from the positive real axis of x. In the example of the soft part of the hard function, singularities appear on the real axis at both positive and negative values. We expect that this would not happen for the full hardfunction. In other words, the singularities present in the resummed soft part of the hard function at positive x must all cancel against similar divergences in the resummed hard part of the hard function. This explains the observations above regarding the radius of convergence of the full hard function versus its soft and hard components, but it applies more generally, also to poles further away from the origin.
To complete the analysis of the radius of convergence in the full hard function we would now like to interpret the numerical values of R obtained in the Padé-based analysis in terms of the colour structures C 1 = (2C A − T 2 t ) and C 2 = (C A − T 2 t ) 9 . We start by recalling that for the soft part of the hard function, R s = 1/C 2 depends on C 2 only. The Padé-based analysis of the 27 colour representation indicates that this is not so for the full hard function. To obtain an analytic expression it proves useful to depart from the actual values of the colour factors corresponding to physically-relevant representations, and simply repeat the Padé approximant analysis for a range of values of C 2 for a fixed C 1 . To this end we plot in figure 13 the numerical values of R emerging from Padé approximants, as a function of C 2 , for fixed values of C 1 (we pick C 1 = 6 and C 1 = −2, corresponding to the singlet and the 27 representation, for easy reference). More precisely, we display in figure  13 the value of 1/R rather than R itself, which makes it easy to recognise the exact linear behaviour. Based on this analysis we deduce the radius of convergence of the full amplitude to be: .  Figure 13. The radius of convergence of the full infrared-renormalized amplitude as a function of the colour operators. In these plots the dots represent the value of 1/R, for the corresponding value of C 1 = (2C A − T 2 t ) and C 2 = (C A − T 2 t ), based on the Padé approximant analysis for = 11 and = 12, as indicated in the plots. We superimpose two linear lines, which determine the dependence of 1/R on the colour operators, as summarised by eq. (6.15).
Returning to the physically-relevant representations, x b end up being closest to the origin (|x b | < |x a |) for the singlet representation where C

Conclusions
In this paper we completed the perturbative calculation of 2 → 2 partonic amplitudes at the next-to-leading logarithmic accuracy in the Regge limit to high loop orders. We focused on the previously-unknown even-signature terms, corresponding to the imaginary part of the amplitude, which vanishes at the leading-logarithmic accuracy. Building upon our previous work in ref. [25] where we determined the infrared singularities, we now computed the finite corrections to the hard amplitude H, which remain after stripping off, or renormalizing, these singularities. We believe that these results -the soft anomalous dimension and hard functions -together exhaust the physical information contained in these partonic amplitudes.
Our results are based on the well-established BFKL evolution equation in momentum space. Since the even amplitude vanishes at the leading-logarithmic order, only the leading-order BFKL evolution kernel was needed in our calculation, and the final formulae apply equally to quark and gluon amplitudes. We exploited the fact, observed in [25], that the two-Reggeon wavefunction is finite. While it is unknown how to diagonalise the BFKL Hamiltonian for arbitrary colour structures beyond the planar limit, we were able to solve the BFKL equation iteratively, treating complementary regions using two different approaches. The first relies on the soft approximation keeping the dimensional regularization parameter finite -the same method we used in ref. [25] to determine the singularities of the amplitude -while the second relies instead on a computation in exactly two transverse dimensions, which captures general hard momentum configurations where both Reggeons carry momenta of the order of the total momentum transfer p 2 = −t. As shown in eqs. (3.9) and (3.13), each separated part of the BFKL-motivated reduced amplitude needs only be calculated to order O( 0 ), and by carefully recombining them we obtained the renormalized amplitude in eq. (5.24). The result passes several consistency checks and agrees with a direct computation in dimensional regularization, which we performed through five loops.
The central new computation in this paper is the iterative solution of the BFKL equation in two dimensions, leading to a simple algorithm to compute the two-Reggeon wavefunction to any order, presented in section 4. The result lives inside a very rigid space of functions: the -loop wavefunction is a linear combination of weight-single-valued harmonic polylogarithms (SVHPLs) of z andz with rational coefficients. The algorithm is formulated as an operation on SVHPLs, and it works by producing differential equations in the holomorphic variable z that can be directly integrated in terms of HPLs of z, to which we subsequently apply the single-value map to recover the actual wavefunction in terms of SVHPLs of z andz.
The hard contribution to the infrared-renormalized amplitude H computed using the two-dimensional method admits a rather complex structure, and its resummation goes beyond the scope of the present paper. This is to be contrasted with the soft contribution, which we could resum to all orders in terms of gamma functions, eq. (3.37), which includes singular as well as finite corrections. The number-theoretical content of the hard contribution is interesting: by construction it is restricted to single-valued multi zeta values (see eqs. (5.19) and (5.21)). The presence of multi zeta values -which make their first appearance at weight 11 involving a single-valued version [27,29] of ζ 5,3,3 -precludes resummation in terms of gamma functions, so the resummed result would clearly be of different nature to that of the soft contribution.
Having obtained explicit analytic expressions for both the two-Reggeon wavefunction and the infrared-renormalized amplitude H to high loop order, it is straightforward to study the results numerically. In section 6 we examine a couple of aspects, first considering the wavefunction and then the infrared-renormalized amplitude. The wavefunction manifests highly regular behaviour as a function of the Reggeon kinematics variables, except for three specific limits. Two of these correspond to the soft limits, z,z → 0 and z,z → ∞, while the third z,z → 1 corresponds to the limit of large internal momentum. The former are described analytically by the soft wavefunction in eq. (3.24), and by definition the hard wavefunction vanishes there, while a peak at large momentum is present in both the soft and hard wavefunctions. Interestingly, there is a significant -but incompletecancellation between these two leading to a more modest peak in the full wavefunction. While this phenomenon does not affect the validity of our results, it would be interesting to independently predict this limit of the wavefunction (extending eq. (6.5)) which could help find simpler numerical approximations. Considering the infrared-renormalized amplitude we focused on one interesting problem, namely the convergence of the perturbative expansion. We find that the O( 0 ) infrared-renormalized amplitude has a finite radius of convergence in the expansion parameter x = Lα s /π. For the soft contribution, where we have a resummed analytic expression, eq. (3.36), this radius of convergence can readily be identified as the first pole of a gamma function, generating asymptotic behaviour ∼ (x(C A − T 2 t )) at high orders, → ∞. The soft contribution is however not physically meaningful on its own, and the complete infrared-renormalized amplitude features a larger radius of convergence, as shown in figure  12 (compare with figures 10 and 11 for the separate soft and hard components). Estimating the convergence radius using Padé approximants for different colour channels, we deduced an empirical formula for the radius of convergence R of the full amplitude in terms of C A and T 2 t , eq. (6.15) above. Interestingly, the pole closest to the origin is always on the negative real axis, leading to an asymptotic behaviour of alternating signs. This matches our physical expectation that the resummed expression should smoothly extrapolate to high energies, corresponding to large positive values of x, and is similar to what was observed previously for non-global logarithms in ref. [51]. It remains for future work to understand the true high-energy (large x = Lα s /π) behaviour.
Let us conclude with a brief summary of the state-of-the-art knowledge of partonic 2 → 2 scattering amplitudes in the Regge limit. With the completion of this work these amplitudes are known in full to NLL accuracy. The signature odd part, corresponding to the exchange of a single Reggeized gluon was already known, and is given by a Regge pole (1.3) with two-loop corrections to the trajectory α (2) g (p 2 ), and suitable impact factors (the former, in particular, was calculated in [52][53][54][55]; it can also be extracted from two-loop calculations of 2 → 2 scattering amplitudes [18]). The signature even part, corresponding to a pair of Reggeized gluons, which generate a Regge cut, was determined here. The next frontier is therefore NNLL accuracy. In the signature-odd sector the first step was taken in ref. [24], where the non-linear Balitsky-JIMWLK equation was used to compute the Regge cut contribution generated through the evolution of three Reggeized gluons and their mixing with one Reggeon through three loops. It is very interesting, and indeedusing the techniques we developed in the present paper -technically feasible, to compute higher-loop corrections in this tower of logarithms. NNLL corrections in the signatureeven sector are in turn simpler and can be deduced from linear BFKL evolution with a NLO kernel [52][53][54][55], supplemented by suitable impact factors. At N 3 LL one expects new phenomena such as the mixing of two and four Reggeon states, which can again be computed using the Balitsky-JIMWLK equation.
Finally, beyond their immediate relevance to the the study of the high-energy limit, the results in this paper can be used to check future multi-loop calculations, and ultimately serve as "boundary data" in a bootstrap programme in which amplitudes are deduced using knowledge of the space of functions, analytic properties, symmetries and special kinematic limits. Such a programme was highly successful in the context of N = 4 supersymmetric Yang-Mills theory, see e.g. [56,57], but also, more recently in the context of the singularity structure of gauge theories including QCD [36]. In both cases, the high-energy limit served as crucial input. the number of nested integrals. The recursion is closed by the weight-1 identities H 0 (z) = log z and H 1 (z) = − log(1 − z).

(A.2)
HPLs form a shuffle algebra and thus obey shuffle product identities where ρ ¡ σ denotes the shuffle of the words ρ and σ.
The indices of a HPL may be shortened by means of a collapsed notation; one replaces strings of zeros followed by a one according to for example H 0,1,0,0,1,1 (z) → H 2,3,1 (z). In the collapsed notation the number of indices is referred to as the depth of the function (while their sum now equals the weight).
Depending on the context it may be useful to view the HPLs as nested sums. One commonly used definition is where we assume the collapsed notation. Note that the aforementioned depth is equal to the number of nested sums. The Taylor series of HPLs, defined by eq. (A.1), whose rightmost index is non-zero, is given by eq. (A.5) with (A.6). Trailing zeros in the indices of a HPL point to logarithmic divergences at z = 0. The log z = H 0 (z) terms can be exposed using the shuffle algebra; one considers H σ (z)H 0 (z) = H σ,0 (z) + . . . + H 0,σ (z) (A.7) and solves for H σ,0 (z). This procedure can be applied recursivly until all trailing zeros are removed. Hence, HPLs can always be written as a series in z and log z. For arguments between 0 and 1 HPLs yield real values. They show branch cuts on the real axis where z ∈ [1, ∞) and are thus multi-valued functions.

B Single-valued harmonic polylogarithms
Single-valued harmonic polylogarithms (SVHPLs) [26] are the class of all branch cut-free, single-valued, combinations of HPLs. Their construction is somewhat involved and we will only provide a short summary here. Further details can be found in e.g. refs. [30][31][32].

(B.2)
For the explicit construction one typically defines two alphabets {x 0 , x 1 } and {y 0 , y 1 } and the corresponding sets of all words X * and Y * formed from the respective alphabet. The letters of the former alphabet directly translate to {0, 1} when they appear as the indices of a (SV)HPL. The letters y 0 , y 1 are related to x 0 , x 1 via Z(y 0 , y 1 )y 1Z (y 0 , y 1 ) −1 = Z(x 0 , where Z is the so-called Drinfeld associator. It is defined as the generating series where with "tilde" and φ defined below eq.(B.5). SVHPLs obey the same shuffle product as HPLs (A.3), namely L ρ (z,z)L σ (z,z) = τ ∈ρ¡σ L τ (z,z). (B.8)

B.1 Holomorphic part and single-value map
SVHPLs are uniquely fixed by their holomorphic part (i.e. their functional dependence on z) and the requirement of single-valuedness. We define the holomorphic part of a function ψ(z,z) as the limit ψ (h) (z) = ψ(z, 0) logz→0 . (B.9) For a given linear combination of SVHPLs taking this limit simply amounts to replacing L σ (z,z) → H σ (z). The dependence onz is reconstructed by the single-value map s ψ (h) (z) = ψ(z,z) (B.10) which is discussed in detail in refs. [27,34]. Again, we restrict ourselves here to stating the (obvious) replacement rule H σ (z) → L σ (z,z) which generates the corresponding singlevalued expression from a linear combination of HPLs of z. As the action of the Hamiltonian H 2d,i (4.11a) removes constant terms from the wavefunction prior to integration we shall not discuss this aspect in the context of eqs. (B.9) and (B.10) here. The interested reader is referred to the above references.
ŝ 1→0+1 to split L σ (z,z) into a sum of 2 n 1 (σ) SVHPLs according to the index rule 1 → 0 + 1. withŝ 0→0+1 defined likeŝ 1→0+1 (B.10) but based on the index rule 0 → 0 + 1. This last step is not strictly necessary but it reduces the amount of data needed to apply these kinds of transformations to a list of SVHPLs at z,z = 1.
Lastly, let us examine the transformation z ↔z and how to related an SVHPL L σ (z, z) to (a sum of) SVHPLs L σ i (z,z). The easy yet computationally heavy way is to translate L σ (z,z) to HPLs, swap z ↔z, extract the holomorphic part by means of eq. (B.9) and finally apply s (B.10). For SVHPLs of weight less or equal to five this might be adequate but at higher weights it becomes inefficient due to the large size of expressions that the translation to HPLs causes. Like in the above examples this step can be avoided altogether.
The procedure relies on knowing the functional dependence of y 1 on the x i , cf. eq. (B.4). Consider the weight-n SVHPL L σ (z,z) with σ = σ 1 , . . . , σ n and swap z ↔z. Then where the "tilde" map was defined below eq. (B.5) and y1(σ) is the coefficient of the product of x 0 and x 1 corresponding to σ, e.g. if σ = 1, 1, 0, 1, 0 then y1(σ) is the coefficient of