---
author:
- 'Simon Caron-Huot,'
- 'Einan Gardi,'
- 'Joscha Reichel,'
- Leonardo Vernazza
bibliography:
- 'main.bib'
title: 'Two-parton scattering amplitudes in the Regge limit to high loop orders'
---

Introduction {#intro}
============

The study of QCD scattering in the Regge limit has been an active area of research for over half a century, e.g. [@Kuraev:1977fs; @Balitsky:1978ic; @Lipatov:1985uk; @Mueller:1993rr; @Mueller:1994jq; @Brower:2006ea; @Moult:2017xpp]. 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\gg-t\gg\Lambda_{\rm QCD}^2$ (see figure \[setup\_fig\]), perturbation theory offers systematic tools to analyse this limit. Central to this is the Balitsky-Fadin-Kuraev-Lipatov (BFKL) evolution equation [@Kuraev:1977fs; @Balitsky:1978ic], 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 [@Balitsky:1995ub; @Balitsky:1998kc; @Kovchegov:1999yj; @JalilianMarian:1996xn; @JalilianMarian:1997gr; @Iancu:2001ad], 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. [@Sotiropoulos:1993rd; @Korchemsky:1993hr; @Korchemskaya:1996je; @Korchemskaya:1994qp; @DelDuca:2001gu; @DelDuca:2013ara; @DelDuca:2014cya; @Bret:2011xm; @DelDuca:2011ae; @Caron-Huot:2013fea; @Caron-Huot:2017fxr; @Caron-Huot:2017zfo]. 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 [@Caron-Huot:2013fea; @Caron-Huot:2017fxr; @Caron-Huot:2017zfo] of $2\to 2$ partonic amplitudes ($qq\to qq$, $gg\to gg$, $qg\to qg$) in QCD and related gauge theories.

![The $t$-channel exchange dominating the high-energy limit, $s\gg -t>0$. The figure also defines our conventions for momenta assignment and Mandelstam invariants. We shall assume that particles 2 and 3 (1 and 4) are of the same type and have the same helicity.[]{data-label="setup_fig"}](./img/setup_fig-crop.pdf)

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 $t$-channel exchange (figure \[setup\_fig\]) of effective degrees of freedom called *Reggeized gluons*. $2\to 2$ amplitudes are conveniently decomposed into *odd* and *even* signature characterising their symmetry properties under $s\leftrightarrow u$ interchange, or crossing symmetry: $$\label{Odd-Even-Amp-Def}
    {\cal M}^{(\pm)}(s,t) = \tfrac12\Big( {\cal M}(s,t) \pm {\cal M}(-s-t,t) \Big)\,,$$ where odd (even) amplitudes ${\cal M}^{(-)}$ (${\cal M}^{(+)}$) are governed by the exchange of an odd (even) number of Reggeized gluons. Furthermore, as shown in ref. [@Caron-Huot:2017fxr], these have respectively *real* and *imaginary* coefficients, when expressed in terms of the natural signature-even combination of logarithms, $$\label{L-def}
    \frac12\left(\log\frac{-s-i0}{-t}+\log\frac{-u-i0}{-t}\right)
    \simeq \log\left|\frac{s}{t}\right| -i\frac{\pi}{2} \equiv L\,.$$

The real part of the amplitude, ${\cal 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 $$\label{Mreal}
    {\cal M}^{(-)}_{\rm LL} = (s/(-t))^{\alpha_g(t)} \times {\cal M}^{\rm tree}$$ where the exponent is the *gluon Regge trajectory* (corresponding to a Regge pole in the complex angular momentum plane), $\alpha_g(t)=\frac{\alpha_s}{\pi} C_A 
\alpha_g^{(1)}(t)+{\cal O}(\alpha_s^2)$, whose leading order coefficient $\alpha_g^{(1)}(t)$ is infrared singular, $\alpha_g^{(1)}(t)\sim \frac{1}{2\epsilon}$ in dimensional regularization with $d=4-2\epsilon$ (see eq. (\[alphag1\]) below). Infrared singularities are well-known to exponentiate, independently of the high-energy limit. Importantly, however, eq. (\[Mreal\]) 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. [@Sotiropoulos:1993rd; @Korchemsky:1993hr; @Korchemskaya:1996je; @Korchemskaya:1994qp; @DelDuca:2001gu; @DelDuca:2013ara; @DelDuca:2014cya; @Bret:2011xm; @DelDuca:2011ae; @Caron-Huot:2013fea; @Caron-Huot:2017fxr; @Caron-Huot:2017zfo].

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 ${\cal O}(\alpha_s^2)$ 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. [@Caron-Huot:2017fxr] 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, ${\cal M}^{(+)}$, extending our work [@Caron-Huot:2017zfo]. 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: $$\label{MevenOneloop}
    {\cal M}^{(+)}_{\rm NLL}\simeq  i\pi 
    \left[\frac{1}{2\epsilon} \frac{\alpha_s}{\pi}
    +{\cal O}\left(\alpha_s^{2} L\right)\right] 
    {\mathbf T}^2_{s-u}  {\cal M}^{\rm tree}\,.$$ Here we suppressed subleading terms in $\epsilon$ as well as multiloop corrections, which take the form $\alpha_s^{\ell} 
L^{\ell-1}$ at $\ell$ 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. (\[MevenOneloop\]) one may observe another salient feature of this tower of corrections, namely the colour structure, which is even under $s\leftrightarrow u$ interchange (${\cal M}^{\rm tree}$ is odd, and so is the operator ${\mathbf T}^2_{s-u}$ acting on it). The first term in the square brackets in (\[MevenOneloop\]) is the exact result in the planar limit; we will be interested in the full series of corrections $\alpha_s^{\ell} L^{\ell-1}$, which are all subleading in the large $N_c$ limit (see the definitions of colour operators in eq. (\[TtTsTu\]) below).

All higher-order corrections, ${\cal O}(\alpha_s^{\ell}L^{\ell-1})$, in (\[MevenOneloop\]) can be described by the well-known ladder graphs, where each additional loop constitute an additional rung in the ladder (see figure \[fig:simplebfkl\] below).

![Sketch of evolution generating ladder graphs in the imaginary part of the amplitude. Considering initially emission from the projectile side only, the 0-loop wavefunction (top left) describes a state involving two reggeized gluons. The Reggeized gluons are both off-shell and are characterized by their transverse momenta $k$ and $p-k$. Each application of the BFKL Hamiltonian (the top row) generates an additional rung in the ladder. Upon integrating the $(\ell-1)$-loop wavefunction with the target one obtains the $\ell$-loop amplitude (bottom row).[]{data-label="fig:simplebfkl"}](./img/simplebfkl-crop.pdf)

Being the leading contributions to the imaginary part of the amplitude, they are particularly important, and clearly at high energies, where $\alpha_s L\sim {\cal O}(1)$, one should aim at an all-order calculation. These corrections, however, do not feature a simple exponentiation pattern as in eq. (\[Mreal\]); 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. [@Caron-Huot:2013fea; @Caron-Huot:2017fxr; @Caron-Huot:2017zfo]. In [@Caron-Huot:2013fea] higher-order terms in eq. (\[MevenOneloop\]) were computed through four loops – the first order where finite contributions appear (see eqs. (28-29) in [@Caron-Huot:2017fxr]). Subsequently, in ref. [@Caron-Huot:2017zfo] 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 $(\ell-1)$ loops, and the corresponding $\ell$-loop $2\to 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\ll p^2 \ll s$. This is the limit used in ref. [@Caron-Huot:2017zfo] 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 \[fig:simplebfkl\]). In the present paper, building upon the computation of the wavefunction in [@Caron-Huot:2017zfo] we introduce a symmetrized solution accounting simultaneously for the two soft limits, $k^2\ll p^2$ and $(p-k)^2\ll  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. [@Brown:2004ugm] and then studied and applied in the context of motivic periods [@Brown:2013gia] and Feynman integrals [@Chavez:2012kn; @Schnetz:2013hqa]. The relevance of this class of functions for gauge-theory amplitudes within the Regge limit [@Pennington:2012zj; @Dixon:2012yy; @DelDuca:2013lma; @Dixon:2014voa; @DelDuca:2016lad; @DelDuca:2018hrv] (and beyond [@Almelid:2017qju; @Dixon:2019lnw]) 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 \[chap:bfkl\] 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 \[soft\] we review the soft approximation developed in [@Caron-Huot:2017zfo] 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 $\alpha_sL$, which resums both infrared-divergent and finite contributions to all loops, within the soft approximation. In section \[2d-bfkl\] 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 \[amplitude\] to compute the corresponding two-to-two scattering amplitude. In section \[numerics\] we perform a numerical study of the resulting wavefunctions and amplitudes, and address the convergence of the perturbative expansion. Finally, in section \[conclusion\] we make some concluding comments and present an outlook for future investigation.

 equation in dimensional regularisation and the $2\to 2$ amplitude {#chap:bfkl}
==================================================================

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 [@Caron-Huot:2013fea; @Caron-Huot:2017fxr]. 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 [@Caron-Huot:2017zfo], 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 \to 2$ scattering amplitudes. We can single out this contribution by introducing a reduced amplitude, in which the one-Reggeon exchange has been removed: \[Mreduced\] \_[ijij]{} e\^[-\_g(t) L [\_t\^2]{}]{} \_[ijij]{}, where $L$ is the signature-even high-energy logarithm defined in eq. (\[L-def\]), ${{\mathbf{T}}_t^2}$ represents the total colour charge exchanged in the $t$ channel (see eq.  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 $$\begin{aligned}
    \alpha_g(t)=\frac{\alpha_s}{\pi} \alpha_g^{(1)}(t)+{\cal O}(\alpha_s^2)\end{aligned}$$ is the *gluon Regge trajectory* introduced already in [eq. ]{}, where the leading-order coefficient in dimensional regularization with $d=4-2\epsilon$ is given by $$\label{alphag1}
    \alpha_g^{(1)}(t)= \frac{B_0}{2\epsilon} 
    \left(\frac{-t}{\mu^2}\right)^{-\epsilon}$$ where $$\label{B0}
    {B_{0}} \equiv {B_{0}}({\epsilon})=e^{{\epsilon}\gamma_{\rm E}}  
    \frac{\Gamma^2(1-{\epsilon})\Gamma(1+{\epsilon}) }{\Gamma(1-2 {\epsilon})}
    = 1 - \frac12 {\epsilon}^2 \zeta_2 -\frac73 {\epsilon}^3\zeta_3 + O({\epsilon}^4)$$ 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. ]{}, thus we focus only on this component in the following. As discussed in [@Caron-Huot:2017zfo], the reduced amplitude takes the form of an integral over the two-Reggeon wavefunction ${\Omega}(p,k)$, as follows: \[ReducedAmpNLL\] [[\_]{}\^[(+)]{}]{}() = -i[\[k\]]{} (p,k) [\_[s-u]{}\^2]{} [\^]{}\_[ijij]{}, where $p^2 = -t$. In [eq. ]{} the integration measure is \[measure\] [\[k\]]{} ( )\^ , and ${{\mathcal{M}}^{\mathrm{(tree)}}}_{ij\to ij}$ represent the tree amplitude, given by [\^]{}\_[ijij]{} = 4\_s (T\_i\^b)\_[a\_1 a\_4]{} (T\_j\^b)\_[a\_2 a\_3]{} \_[\_1\_4]{}\_[\_2\_3]{}, where $\lambda_i$ for $i=1$ through $4$ are helicity indices. The colour operator ${{\mathbf{T}}_{s-u}^2}$ in eq.  acts on ${{\mathcal{M}}^{\mathrm{(tree)}}}_{ij\to ij}$ and it is defined in terms of the usual basis of quadratic Casimirs corresponding to colour flow through the three channels [@Dokshitzer:2005ig; @DelDuca:2011ae]: \[TtTsTu\] [\_[s-u]{}\^2]{} [with]{} {

[c]{} \_s = \_1+\_2=-\_3-\_4,\
\_u = \_1+\_3=-\_2-\_4,\
\_t = \_1+\_4=-\_2-\_3,

. where ${\mathbf{T}}_i$ is the colour-charge operator [@Catani:1998bh] associated with parton $i$.

The BFKL equation [@Kuraev:1977fs; @Balitsky:1978ic] for the wavefunction ${\Omega}(p,k)$ in eq. (\[ReducedAmpNLL\]) takes the form $$\label{BFKL_evolution}
    \frac{d}{dL}\Omega(p,k)=
    \frac{\alpha_s B_0(\epsilon)}{\pi} \hat{H} \Omega(p,k)\,,$$ where $L$ is the high-energy logarithm (\[L-def\]) and where the Hamiltonian takes the form [@Caron-Huot:2017zfo] \[Hdef1\] = [(2[C\_A]{}-[\_t\^2]{})]{} [\_]{}+ [([C\_A]{}-[\_t\^2]{})]{} [\_]{}, where two independent colour factors come along with two different operations:

\[Hamil\] $$\begin{aligned}
        {{\hat{H}}_{\mathrm{i}}}\, \Psi(p,k) &= \int {[{\mathrm{D}}k']}\, f(p,k,k') \left[ \Psi(p,k') - \Psi(p,k) \right]  \label{Him}, \\
        {{\hat{H}}_{\mathrm{m}}}\, \Psi(p,k) &= J(p,k) \, \Psi(p,k)\,   \label{eq:Hm}\,.
    \end{aligned}$$

The function $f(p,k,k')$ in eq. (\[Him\]) represents the evolution kernel \[bfkl-kernel\] f(p,k,k’) + - , and $J(p,k)$ in eq. (\[eq:Hm\]) is defined by $$\begin{aligned}
  \label{Jp-def2}
    J(p,k) = \frac{1}{2{\epsilon}} + \int {[{\mathrm{D}}k']}\, f(p,k,k')
    = \frac{1}{2{\epsilon}} \left[2- {\left(\frac{p^2}{k^2}\right)}^{{\epsilon}} - {\left(\frac{p^2}{(p-k)^2}\right)}^{{\epsilon}} \right].\end{aligned}$$ While it is unknown how to diagonalise this $d$-dimensional Hamiltonian, we may invoke a perturbative solution [@Caron-Huot:2013fea; @Caron-Huot:2017zfo] by expanding the wavefunction in the strong coupling constant: \[OmegaEven\] (p,k) &=& \_[=1]{}\^ ( )\^ L\^[-1]{} [\^[([-1]{})]{}]{}(p,k), where we set the renormalisation scale equal to the momentum transfer, $\mu^2 = -t = p^2$.

Substituting the expanded form of the wavefunction in (\[OmegaEven\]) into the BFKL evolution equation (\[BFKL\_evolution\]) one deduces that \[Hdef0\] [\^[([-1]{})]{}]{}(p,k) = [\^[([-2]{})]{}]{}(p,k), where ${\hat{H}}$ is the BFKL hamiltonian of eq. (\[Hdef1\]), that is, the wavefunction at any given order is found by repeated application of the Hamiltonian, where the initial condition in our normalization is simply \[0th-wavefunction-tilde\] [\^[([0]{})]{}]{}(p,k) = 1.

Next, let us consider the on-shell $2\to 2$ amplitude. Substituting the expanded wavefunction (\[OmegaEven\]) into (\[ReducedAmpNLL\]) we readily obtain the following expansion \[MhatEven\] [[\_]{}\^[(+)]{}]{}( ) &=& \_[=1]{}\^( )\^ L\^[-1]{} [[\_]{}\^[(+,)]{}]{}, with \[ReducedAmpNLL2\] [[\_]{}\^[(+,)]{}]{} = -i [\[k\]]{} [\^[([-1]{})]{}]{}(p,k) [\_[s-u]{}\^2]{}[\^]{} . Namely, integrating over the $(\ell-1)$-th order contribution to the wavefunction yields the $\ell$-th order contribution to the amplitude. A graphical illustration of eq.  is provided in figure \[fig:bfklwfamp\]. As discussed in the introduction, because of evolution, the amplitude at accuracy can be represented as a ladder. At order $\ell$ it is obtained by closing the ladder and integrating the wavefunction of order $(\ell-1)$ over the resulting loop momentum, according to eq. . The wavefunction ${{\Omega}^{({\ell-1})}}(p,k)$ in turn is obtained by applying once the leading-order evolution kernel to the wavefunction of order $(\ell-2)$. Graphically, this operation corresponds to adding one rung to the ladder.

![Graphical representation of the amplitude at accuracy, as obtained through evolution. The addition of one rung corresponds to applying once the leading-order evolution on the wavefunction of order $(\ell-2)$. This gives the wavefunction at order $(\ell-1)$, according to eq. . Closing the ladder and integrating over the resulting loop momentum gives the reduced amplitude, according to eq. .[]{data-label="fig:bfklwfamp"}](./img/bfklwfamp-crop.pdf)

Inspecting [eqs.  and ]{} 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

$$\begin{aligned}
        \label{J_im_general}
        {\Omega}_{\mathrm{i},w}(p,k) &\equiv \int {[{\mathrm{D}}k']}f(p,k,k') \left[ {\Omega}_{w}(p,k') - {\Omega}_{w}(p,k) \right],  \\
        {\Omega}_{\mathrm{m},w}(p,k) &\equiv J(p,k) \, {\Omega}_{w}(p,k) \label{eq:Wmdef}, 
    \end{aligned}$$

where ${\Omega}_{\varnothing}(p,k) \equiv 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.  and (\[eq:Hm\]), respectively. In this notation the first four orders of the wavefunction read, for instance, $$\begin{aligned}
    {{\Omega}^{({1})}}(p,k) &= {({C_A}-{{\mathbf{T}}_t^2})}{\Omega}_{\rm m},  \label{WavefunctionTwoLoops-b} \\
    {{\Omega}^{({2})}}(p,k) &= {({C_A}-{{\mathbf{T}}_t^2})}^2 {\Omega}_{\rm m,m} + {(2{C_A}-{{\mathbf{T}}_t^2})}{({C_A}-{{\mathbf{T}}_t^2})}{\Omega}_{\rm i,m},  \label{WavefunctionTwoLoops-b2} \\ \nn
    {{\Omega}^{({3})}}(p,k) &= {({C_A}-{{\mathbf{T}}_t^2})}^3 {\Omega}_{\rm m,m,m} + {(2{C_A}-{{\mathbf{T}}_t^2})}{({C_A}-{{\mathbf{T}}_t^2})}^2 \left( {\Omega}_{\rm i,m,m} + {\Omega}_{\rm m,i,m} \right) \\
    & +\,  {(2{C_A}-{{\mathbf{T}}_t^2})}^2 {({C_A}-{{\mathbf{T}}_t^2})}{\Omega}_{\rm i,i,m},  \label{eq:wf3loops}  \\ \nn
    {{\Omega}^{({4})}}(p,k) &= {({C_A}-{{\mathbf{T}}_t^2})}^4 {\Omega}_{\rm m,m,m,m} \\ \nn
    & +\, {(2{C_A}-{{\mathbf{T}}_t^2})}{({C_A}-{{\mathbf{T}}_t^2})}^3 \left( {\Omega}_{\rm m,m,i,m} + {\Omega}_{\rm m,i,m,m}  + {\Omega}_{\rm i,m,m,m}  \right) \\ \nn
    & +\,  {(2{C_A}-{{\mathbf{T}}_t^2})}^2 {({C_A}-{{\mathbf{T}}_t^2})}^2 \left( {\Omega}_{\rm m,i,i,m} + {\Omega}_{\rm i,m,i,m} + {\Omega}_{\rm i,i,m,m} \right) \\
    & +\,  {(2{C_A}-{{\mathbf{T}}_t^2})}^3 {({C_A}-{{\mathbf{T}}_t^2})}{\Omega}_{\rm i,i,i,m}. \label{eq:wf4loops}\end{aligned}$$

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 evolution in figure \[fig:bfklwfamp\]. This implies \[left-right-symmetry\] [\^[()]{}]{}(p,k) = [\^[()]{}]{}(p,p-k), which can be easily verified by showing that the functions $f(p,k,k')$ in , $J(p,k)$ in and ${{\Omega}^{({0})}}(p,k)$ in obey the same symmetry. This symmetry property will become handy in section \[soft\], making it possible to capture simultaneously both soft limits, $k^2\to 0$ and $(p-k)^2\to 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{\epsilon}$ 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 $$\begin{aligned}
    \label{triang1} \nn
    {\Omega}_{\rm i,m}  &\ni \int {[{\mathrm{D}}k']}\, 
    \frac{(p-k)^2}{(p-k')^2 (k-k')^2} {\left(\frac{p^2}{(k')^2}\right)}^{{\epsilon}}, \\ 
    {\Omega}_{\rm i,i,m}  &\ni  \int {[{\mathrm{D}}k']}[Dk''] \, 
    \frac{k^2 (p-k'')^2}{(k'')^2  (p-k')^2 (k-k'')^2 (k'-k'')^2} {\left(\frac{p^2}{(k')^2}\right)}^{{\epsilon}},\end{aligned}$$ which are represented respectively in figure \[fig:triang\] (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 ${\epsilon}$ in eq. .[]{data-label="fig:triang"}](./img/Wavefunction-loops.pdf){width="72.00000%"}

The amplitude at order $\ell$ is obtained upon integrating the wavefunction of order $\ell-1$, as indicated in [eq. ]{}. 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 [@Caron-Huot:2017zfo] \[JiJmSymAllOrders\] [\[k\]]{} [\_]{}\_[w]{}(p,k) = [\[k\]]{} \_[,w]{}(p,k) = 0. Furthermore, in the notation of eqs.  and reversal of the rungs directly translates to the reversal of the indices of the wavefunction. The target-projectile symmetry thus guarantees the equality [\[k\]]{} \_[a\_1,…,a\_n]{}(p,k) = [\[k\]]{} \_[a\_n,…,a\_1]{}(p,k). 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 \[fig:5loopex\] (a). Up to four loops one obtains [@Caron-Huot:2017zfo] $$\begin{aligned}
    \label{ReducedAmpNLL2-one-loop}
    {{{\hat{{\mathcal{M}}}}_{\mathrm{NLL}}}^{(+,{1})}} &= i\pi \frac{{B_{0}}}{2{\epsilon}} {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}},  \\
    \label{ReducedAmpNLL2-two-loop} \nn
    {{{\hat{{\mathcal{M}}}}_{\mathrm{NLL}}}^{(+,{2})}} &= i\pi \frac{({B_{0}})^2}{2} \left[ \frac{1}{(2{\epsilon})^2} 
    + \frac{9\zeta_3}{2}{\epsilon}+ \frac{27\zeta_4}{4}{\epsilon}^2 
    + \frac{63\zeta_5}{2}{\epsilon}^3 + {\mathcal{O}}({\epsilon}^4) \right]  \\
    &\hspace{55mm} \times {({C_A}-{{\mathbf{T}}_t^2})}{{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}, \\ \nn
    \label{ReducedAmpNLL2-three-loops}
    {{{\hat{{\mathcal{M}}}}_{\mathrm{NLL}}}^{(+,{3})}} &= i\pi \frac{{B_{0}}^3}{3!} \left[ \frac{1}{(2{\epsilon})^3} 
    - \frac{11\zeta_3}{4} - \frac{33\zeta_4}{8}{\epsilon}- \frac{357\zeta_5}{4}{\epsilon}^2 + {\mathcal{O}}({\epsilon}^3) \right] \\
    & \hspace{55mm} \times {({C_A}-{{\mathbf{T}}_t^2})}^2 {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}, \\ \nn
    \label{ReducedAmpNLL2-four-loops}
    {{{\hat{{\mathcal{M}}}}_{\mathrm{NLL}}}^{(+,{4})}} &= i\pi \frac{{B_{0}}^4}{4!} 
    \bigg\{ {({C_A}-{{\mathbf{T}}_t^2})}^3 \left( \frac{1}{(2{\epsilon})^4} 
    + \frac{175\zeta_5}{2}{\epsilon}+ {\mathcal{O}}({\epsilon}^2) \right) \\ 
    & \hspace{8mm} + {C_A}{({C_A}-{{\mathbf{T}}_t^2})}^2 \left( -\frac{\zeta_3}{8{\epsilon}} 
    - \frac{3}{16}\zeta_4 - \frac{167\zeta_5}{8}{\epsilon}+ {\mathcal{O}}({\epsilon}^2) \right) \bigg\} {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}. \end{aligned}$$ A thorough discussion of the target-projectile symmetry, and its effect on the colour structure of the amplitude has been given in [@Caron-Huot:2017zfo], 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. [@Caron-Huot:2013fea; @Caron-Huot:2017zfo] to higher loops is provided by numerical integration methods. In particular, we find sector decomposition as implemented in `pySecDec`/`SecDec` [@Carter:2010hi; @Borowka:2017esm] 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 $\ell$-loop amplitude: Assigning $o({\epsilon}) = -1$, $o(\pi) = 1$ and $o(\zeta_n) = n$ one sees that the terms of the $\ell$-loop amplitude are uniformly of weight $o({{{\hat{{\mathcal{M}}}}_{\mathrm{NLL}}}^{(+,{\ell})}} ) = \ell$. We can hence deduce which zeta numbers (or powers of $\pi$) may appear at any given order in ${\epsilon}$.

Another observation facilitates this procedure at five loops; after dividing the $\ell$-loop amplitude by $B_0^\ell$ there are no occurrences of $\zeta_2 = \pi^2/6$ up to four loops, see e.g. the ${\mathcal{O}}({\epsilon})$ terms of eq. . If we assume this absence of $\zeta_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, $\zeta_5$, whereas $\zeta_3 \zeta_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 \[sec:finiteamp\].

To obtain the five-loop amplitude ${{{\hat{{\mathcal{M}}}}_{\mathrm{NLL}}}^{(+,{5})}}$ we integrate the four-loop wavefunction ${{\Omega}^{({4})}}(p,k)$ of  according to eq. . 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 \[fig:5loopex\]. 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: \[eq:5loopex\] \~ .

![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. .[]{data-label="fig:5loopex"}](./img/Amplitude-loops.pdf){width="76.00000%"}

After combining all contributions (and reconstructing the zeta numbers in case of the numerical results) we find $$\begin{gathered}
    {{{\hat{{\mathcal{M}}}}_{\mathrm{NLL}}}^{(+,{5})}} = i\pi \frac{{B_{0}}^5}{5!} 
    \left\{ {({C_A}-{{\mathbf{T}}_t^2})}^4 \left( \frac{1}{32 {\epsilon}^5} - \frac{53 \zeta_5}{2} \right) \right. \\
    \left. + {C_A}{({C_A}-{{\mathbf{T}}_t^2})}^3 \left( -\frac{\zeta_3}{16 {\epsilon}^2} - \frac{3 \zeta_4}{32 {\epsilon}} 
    + \frac{253 \zeta_5}{16} \right) - \frac{5}{2} {C_A}^2 {({C_A}-{{\mathbf{T}}_t^2})}^2 \zeta_5 \right\} 
    {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}. \label{eq:m5num}\end{gathered}$$ This result will serve as a consistency check for our computation below.

The soft approximation {#soft}
======================

In section \[chap:bfkl\] we have shown how the two-Reggeon contribution to the two-parton scattering amplitude is conveniently described in terms of the reduced amplitude $\hat {\cal M}$. The latter is defined in [eq. ]{} by (multiplicatively) removing the single-Reggeon effect from the full amplitude ${\cal M}$. This allowed us to use BFKL evolution to express the two-Reggeon contribution to $\hat {\cal M}$ in terms of iterated integrals. Beyond four loops these integrals become difficult to evaluate exactly in $d = 4-2{\epsilon}$ 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 ${\cal M}$ can be factorized, ${\cal M}={\bf Z} 
{\mathcal{H}}$, where the “infrared renormalization” factor ${\bf Z}$ captures all divergences (which famously exponentiate in terms of the soft anomalous dimension, see e.g. [@Sterman:1995fz; @Collins:1989gx; @Korchemskaya:1994qp; @Catani:1998bh; @Aybat:2006mz; @Sterman:2002qn; @Gardi:2009qi; @Becher:2009cu; @Becher:2009qa; @Almelid:2015jia; @Almelid:2017qju]) while the infrared-renormalized amplitude ${\mathcal{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 ${\cal 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* ${\mathcal{H}}$.

In Ref. [@Caron-Huot:2017zfo] we have shown that the soft anomalous dimension associated with the signature-even amplitude, or indeed the relevant infrared renormalization factor ${\bf Z}$, can be computed to all orders by evaluating the reduced amplitude $\hat{\mathcal{M}}$ to ${\mathcal{O}}({\epsilon}^{-1})$. Similarly, we are going to show now (section \[sec:IRfact\]) that the infrared-renormalized amplitude ${\mathcal{H}}$ (in four dimensions) can be completely determined from the reduced amplitude $\hat{\mathcal{M}}$, evaluated at the same accuracy, i.e. to ${\mathcal{O}}({\epsilon}^0)$. This, along with the fact that the corresponding wavefunction $\Omega$ 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 \[sec:soft\_hard\_split\] 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 \[2d-bfkl\]). These two wavefunction components will subsequently serve to compute the corresponding soft and hard contributions to the reduced amplitude $\hat{\mathcal{M}}$ to the required order, ${\mathcal{O}}({\epsilon}^0)$. In section \[sec:softwave\] we review the main results of Ref. [@Caron-Huot:2017zfo] 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 \[sec:softAmpl\] to compute the corresponding ${\mathcal{O}}({\epsilon}^0)$ contributions to the reduced amplitude. Finally, in section \[hardFdef\] we make use of the results of sections \[sec:IRfact\] and \[sec:softAmpl\] to evaluate the ${\mathcal{O}}({\epsilon}^0)$ soft contributions to the infrared-renormalized amplitude ${\mathcal{H}}$.

Infrared factorisation in the high-energy limit\[sec:IRfact\]
-------------------------------------------------------------

According to the infrared factorisation theorem (see e.g. [@Sterman:1995fz; @Collins:1989gx; @Korchemskaya:1994qp; @Catani:1998bh; @Aybat:2006mz; @Sterman:2002qn; @Gardi:2009qi; @Becher:2009cu; @Becher:2009qa; @Almelid:2015jia; @Almelid:2017qju]), infrared singularities of an amplitude ${\mathcal{M}}$ are multiplicatively renormalised by a factor ${\bf Z}$, \[IRfacteq\] ({p\_i},, [\_s]{}() ) = [**Z**]{} ({p\_i},, [\_s]{}() ) ({p\_i},, [\_s]{}() ), such that the infrared-renormalized amplitude ${\cal H}$ is finite as ${\epsilon}\to 0$. We use a minimal subtraction scheme, where the renormalisation factor ${\bf Z}$ consists of pure poles. It is then given explicitly as the path-ordered exponential of the soft anomalous dimension: \[RGsol\] [**Z**]{} ({p\_i},, [\_s]{}() ) = [P]{} { -\_0\^ [****]{} ({p\_i},, [\_s]{}() ) }, where, to the accuracy needed in this paper, we can restrict to tree-level running coupling: ${\alpha_s}(\lambda) = {\alpha_s}(p)
\left(p^2/\lambda^2\right)^{{\epsilon}}$. Given that ${\bf Z}$ was determined in Ref. [@Caron-Huot:2017zfo] to NLL accuracy in the high-energy logarithm, our goal here is to determine the infrared-renormalized amplitude ${\cal H}$ to the same accuracy. Thus we need to specialise [eq. ]{} to the high-energy limit. Recalling that in this limit the amplitude splits naturally into even and odd components under the $s \leftrightarrow u$ signature symmetry, we may focus directly on the even component (the odd component was analysed already in [@Caron-Huot:2017fxr]): \[MtoHeven\] \^[(+)]{}\_[NLL]{} = [**Z**]{}\^[(-)]{}\_[NLL]{} \^[(-)]{}\_[LL]{} + [**Z**]{}\^[(+)]{}\_[LL]{} \^[(+)]{}\_[NLL]{}. Our final goal is to determine ${\mathcal{H}}^{(+)}_{\rm NLL}$. Let us begin by inverting (\[MtoHeven\]), i.e. \[getH\] \^[(+)]{}\_[NLL]{} = - ([**Z**]{}\^[-1]{})\^[(+)]{}\_[LL]{} [**Z**]{}\^[(-)]{}\_[NLL]{} \^[(-)]{}\_[LL]{} +([**Z**]{}\^[-1]{})\^[(+)]{}\_[LL]{} \^[(+)]{}\_[NLL]{}. In [eq. ]{} both the leading- and next-to-leading logarithmic renormalisation factors are known: ${\bf Z}^{(+)}_{\rm LL}$, and hence also $\big({\bf Z}^{-1}\big)^{(+)}_{\rm LL}$ is easily determined from the single-Reggeon exchange, see [eqs.  and ]{}: \[Zll\] [**Z**]{}\^[(+)]{}\_[LL]{} = e\^[ \_t\^2]{} ([**Z**]{}\^[-1]{})\^[(+)]{}\_[LL]{}= e\^[- \_t\^2]{}, where we defined $x\equiv \frac{\alpha_s}{\pi} L$. The factor ${\bf Z}^{(-)}_{\rm NLL}$ was determined to all orders in perturbation theory in [@Caron-Huot:2017zfo]: comparing eqs. (4.12), (4.14) and (4.17) there we express ${\bf Z}^{(-)}_{\rm NLL} $ as \[Znll\] e\^[- \_t\^2]{} [**Z**]{}\^[(-)]{}\_[NLL]{} = i ( 1 - R() )\^[-1]{} [\_[s-u]{}\^2]{}|\_[poles]{}, where the function $R({\epsilon})$ reads R() &=& -1\
\[Rdef\] &=& -2\_3 \^3 -3\_4 \^4 -6\_5 \^5 -(10 \_6-2\^2\_3 ) \^6 + [O]{}(\^7), see also eq. (3.16) of [@Caron-Huot:2017zfo]. In [eq. ]{} the factor $\exp[- (x {\mathbf{T}}_t^2)/(2{\epsilon})]$ on the l.h.s. is left there, because it corresponds directly to the factor $\big({\bf Z}^{-1}\big)^{(+)}_{\rm LL}$ appearing in [eq. ]{} to the left of ${\bf Z}^{(-)}_{\rm NLL}$. Notice also that the $-1$ in the numerator of the first fraction on the r.h.s. of [eq. ]{} can actually be removed, given that we need to consider only the poles originating from [eq. ]{}, and this term contributes only at ${\mathcal{O}}({\epsilon}^0)$, given that the second $\epsilon$-dependent factor is regular, i.e. $\big[ 1 - C_A/(C_A -{{\mathbf{T}}_t^2}) 
\, R({\epsilon}) \big]^{-1} = 1 + {\mathcal{O}}({\epsilon}^3)$.

[Eq. ]{} contains also the leading-logarithmic infrared-renormalized amplitude ${\mathcal{H}}^{(-)}_{\rm LL}$, which, as in case of ${\bf Z}^{(+)}_{\rm LL}$, is determined by single-Reggeon exchange, compare again with [eqs.  and ]{}: \[Hll\] \^[(-)]{}\_[LL]{} = e\^[ x C\_A]{} \^[(tree)]{}, where we have substituted ${{\mathbf{T}}_t^2}\to C_A$, given that in ${\mathcal{H}}^{(-)}_{\rm LL}$ the operator ${{\mathbf{T}}_t^2}$ acts on the tree-level amplitude.

By this point we collected all ingredients needed to explicitly write down the first term in [eq. ]{}. The only missing term on the r.h.s. of this equation is thus the even amplitude itself, ${\mathcal{M}}^{(+)}_{\rm NLL}$. As explained above, in order to determine ${\mathcal{M}}^{(+)}_{\rm NLL}$ by means of BFKL evolution, we wish to express it in terms of the reduced amplitude ${\cal \hat M}^{(+)}_{\rm NLL}$ of [eq. ]{}. Substituting eqs. (\[Mreduced\]), (\[Zll\]) and (\[Hll\]) into [eq. ]{} we get $$\label{getH2}
    {\mathcal{H}}^{(+)}_{\rm NLL} 
    = - e^{-\frac{x}{2{\epsilon}} {\mathbf{T}}_t^2} {\bf Z}^{(-)}_{\rm NLL} 
    \, e^{\frac{{B_{0}}({\epsilon})-1}{2{\epsilon}} \, x \, C_A} {\mathcal{M}}^{\rm (tree)}
    +e^{\frac{B_0({\epsilon})-1}{2{\epsilon}} \, x \, {\mathbf{T}}_t^2} {\cal \hat M}^{(+)}_{\rm NLL}\,\,,$$ where the factor $e^{-\frac{x}{2{\epsilon}} {\mathbf{T}}_t^2} {\bf Z}^{(-)}_{\rm NLL}$ of (\[Znll\]) can be readily substituted as well (this will be done in section \[hardFdef\]). Eq. (\[getH2\]) is an important step because (given that $B_0({\epsilon})-1 = {\cal O}({\epsilon}^2)$, eq. (\[B0\])) it clearly shows that the hard function ${\mathcal{H}}^{(+)}_{\rm NLL} $ at $\epsilon\to 0$ is completely determined once the BFKL-motivated reduced amplitude ${\cal \hat M}^{(+)}_{\rm NLL}$ is known to ${\mathcal{O}}({\epsilon}^0)$, which is the result anticipated at the beginning of this section. With this in mind, we proceed to compute ${\cal \hat M}^{(+)}_{\rm NLL}$ to ${\mathcal{O}}({\epsilon}^0)$.

Soft and hard wavefunction and amplitude {#sec:soft_hard_split}
----------------------------------------

Our strategy to compute the finite part of the reduced amplitude ${\cal \hat M}^{(+)}_{\rm NLL}$ at higher orders is to separate soft and hard components of the wavefunction and truncate the latter to two transverse dimensions (${\epsilon}=0$), where BFKL evolution is much more tractable (see section \[2d-bfkl\]).

As demonstrated in ref. [@Caron-Huot:2017zfo], the soft limit of the wavefunction, where one of the two Reggeons has a small momentum, e.g., $k^2\ll (p-k)^2\simeq p^2$, fully determines all the singular parts in $\epsilon$. This was used to obtain the all-order result for the renormalisation factor ${\bf Z}^{(-)}_{\rm NLL} $ in [eq. ]{}. In addition, the soft limit generates some ${\cal O}({\epsilon}^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 ${\cal O}({\epsilon}^0)$ terms a clear separation between the two regions is necessary. We choose to do this at the level of the wavefunction ${\Omega}(p,q)$. Recall that ${\Omega}(p,q)$ is a finite function[^2] of ${\epsilon}$ [@Caron-Huot:2017zfo], i.e. any singularities in the reduced amplitude are generated through the final integration over the wavefunction in (\[ReducedAmpNLL\]). To proceed we split the wavefunction into two terms: \[OmegaSplitDef\] (p,k) = [\_]{}(p,k) + [\_]{}(p,k), such that the second term, the hard component, vanishes in soft limits: \[h\_definition\] \_[k0]{}[\_]{}(p,k) = \_[kp]{}[\_]{}(p,k) =0. It then follows from (\[ReducedAmpNLL\]) that no singularities can be generated upon integrating ${{\Omega}_{\mathrm{h}}}(p,k)$ (i.e. all the singularities in ${{{\hat{{\mathcal{M}}}}_{\mathrm{NLL}}}^{(+)}}$ are generated upon integrating ${{\Omega}_{\mathrm{s}}}(p,k)$) and hence only the ${\epsilon}\to 0$ limit of ${{\Omega}_{\mathrm{h}}}$ contributes to the finite part of the reduced amplitude. Denoting the wavefunction in this limit as \[WhardTwod\] [\_]{}\^[([2d]{})]{}(p,k) \_[0]{} [\_]{} = [\_]{}(p,k) - [\_]{}\^[([2d]{})]{}(p,k), the reduced amplitude (\[ReducedAmpNLL\]), through order ${\cal O}({\epsilon}^0)$, is then given as a sum of soft and hard components: \[eq:redampSplit\] [[\_]{}\^[(+)]{}]{}() = \^[(+)]{}\_[NLL,s]{}() + \^[(+)]{}\_[NLL,h]{}() with

\[Msh\] $$\begin{aligned}
    \hat {\cal M}^{(+)}_{\rm NLL,s}\left(\frac{s}{-t}\right) 
    &=-i\pi  \int {[{\mathrm{D}}k]}\,\frac{p^2}{k^2(p-k)^2} {{\Omega}_{\mathrm{s}}}(p,k)\,\,  {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}_{ij\to ij}\,, \\
\label{Mhard}
    \hat {\cal M}^{(+)}_{\rm NLL,h}\left(\frac{s}{-t}\right) 
    &= -i\pi\lim_{{\epsilon}\to 0}\int {[{\mathrm{D}}k]}\,\frac{p^2}{k^2(p-k)^2}\,
    {{\Omega}_{\mathrm{h}}}^{({\rm 2d})}(p,k)\,\, {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}_{ij\to ij}\, .\end{aligned}$$

Equations (\[eq:redampSplit\]) and (\[Msh\]) 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 ${\epsilon}$, but we are allowed to truncate the hard wavefunction to ${\cal O}({\epsilon}^0)$. Note that in (\[Mhard\]) we have already substituted the two-dimensional limit of the hard wavefunction, so taking the ${\epsilon}\to 0$ limit simply amounts to taking the integration momentum $k$ to be two-dimensional. These finite integrals will be done in section \[amplitude\].

Let us briefly summarise our plan for the reminder of this section. After reviewing the main arguments of [@Caron-Huot:2017zfo], our aim in section \[sec:softwave\] is to present a symmetrized version of the soft wavefunction in dimensional regularization, eq. (\[Well-1-ansatz-sym\]), which simultaneously captures the two regions where either of the two Reggeons is soft. We then extract the ${\cal O}(\epsilon^0)$ terms in the wavefunction and resum them; these will be used in section \[amplitude\] to determine the two-dimensional hard wavefunction ${{\Omega}_{\mathrm{h}}}^{({\rm 2d})}(p,k)$ from the full one according to eq. (\[WhardTwod\]). Subsequently in section \[sec:softAmpl\] we use the soft wavefunction, computed to all-orders in $\epsilon$, to determine the corresponding contributions to the reduced $2\to 2$ amplitude. We also present an analytic formula resumming these corrections in eq. (\[ReducedAmpNLLresum2\]). Finally, in section \[hardFdef\] we determine the soft wavefunction contribution to the infrared-renormalized amplitude ${\mathcal{H}}^{(+)}_{\rm NLL} $ using eq. (\[getH2\]).

The soft wavefunction {#sec:softwave}
---------------------

The central property of the wavefunction $\Omega(p,k)$ highlighted in [@Caron-Huot:2017zfo] and already mentioned above, is the fact that it is finite for ${\epsilon}\to 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 (\[ReducedAmpNLL\]), and originate from the soft limits $k \to 0$ and $k \to p$ of $\Omega(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 ${{\Omega}^{(j)}}(p,k)$, with $k$ soft, implies that the momentum $k'$ in ${{\Omega}^{(j-1)}}(p,k')$, which has one rung fewer, can also be taken soft, $k' \to 0$, without affecting the result for ${{\Omega}^{(j)}}(p,k)$. In other words, starting with ${{\Omega}^{(j)}}(p,k)$ where $k$ is soft is equivalent to considering the entire side rail of the ladder consisting of soft momenta, $k'$, $k''$, $\ldots\to 0$. Similarly, starting with $k \to p$ implies that all momenta $(p-k)$, $(p-k')$, $(p-k'')$, $\ldots$, are soft. The symmetry of [eq. ]{}, then, implies that $\Omega(p,k)$ in the two limits $k \to 0$ and $k \to p$ must be the same.

In the soft limit the BFKL hamiltonian becomes [@Caron-Huot:2017zfo] \^[(-1)]{}\_s(p,k)&=&H\_s \^[(-2)]{}\_s(p,k),\
H\_s (p,k) &=& (2C\_A-[\_t\^2]{})\
&& + (C\_A-[\_t\^2]{}) J\_s(p,k) (p,k), \[softH\] where \[JpSoft\] J\_s(p,k) = , is the soft approximation of [eq. ]{}. One finds that the wavefunction becomes a polynomial in $ \xi \equiv (p^2/k^2)^{{\epsilon}}$, i.e., the soft limit turns BFKL evolution into a one-scale problem. The integrals involved in [eq. ]{} are simple bubble integrals of the type \[bubbleGeneral1\] [()]{}\^[n ]{} = -[()]{}\^[(n+1)]{}, where the integration measure is given in eq. (\[measure\]), and the class of bubble functions ${B_{n}}({\epsilon})$ is \[bubbleGeneral2\] [B\_[n]{}]{}() = e\^[\_[E]{}]{} . Note that $B_0$ of (\[B0\]) appearing in the gluon Regge trajectory and in the measure (\[measure\]) corresponds to the special case of (\[bubbleGeneral2\]) with $n=0$.

Using eq. (\[bubbleGeneral1\]) one can write the action of the soft Hamiltonian (\[softH\]) on any monomial ($m\geq 0$): \[softHpower\] H\_s \^m &=& ((1-)(C\_A-[\_t\^2]{}) + [B\_[m]{}]{}()(2C\_A-[\_t\^2]{}))\
&=& (\^m - \^[m+1]{} ),where we have introduced the notation \[bubblehat\] B\_n() 1- = 2 n (2 + n) \_3 \^3 + 3 n (2 + n) \_4 \^4 +…. By making repeated use of [eq. ]{} one finds that the wavefunction at order $(\ell-1)$ can be expressed in a closed-form, as follows [@Caron-Huot:2017zfo]: \[Well-1-ansatz\] [\_s\^[(-1)]{}]{}(p,k) = \_[n=0]{}\^[-1]{} (-1)\^n [()]{}\^[n]{} \_[m=0]{}\^[n-1]{} {1 - [B\_[m]{}]{}() }. As discussed in [@Caron-Huot:2017zfo], 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. ]{} is perfectly valid in the soft limit, it breaks explicitly the symmetry of [eq. ]{} 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. ]{} under $k \leftrightarrow (p-k)$: $$\begin{aligned}
 \label{Well-1-ansatz-sym} \nn
    {{\Omega}_s^{(\ell-1)}}(p,k) &=
    \frac{{({C_A}-{{\mathbf{T}}_t^2})}^{\ell-1}}{(2{\epsilon})^{\ell-1}} \sum_{n=0}^{\ell-1} 
    (-1)^n \binom{\ell-1}{n} {\left(\frac{p^2}{k^2}\right)}^{n{\epsilon}} 
    \bigg(\frac{p^2}{(p-k)^2}\bigg)^{n{\epsilon}} 
    \\ &\hspace{3.0cm} \times
    \, \prod_{m=0}^{n-1} 
    \left\{1 - {\hat B_{m}}({\epsilon}) \frac{{(2{C_A}-{{\mathbf{T}}_t^2})}}{{({C_A}-{{\mathbf{T}}_t^2})}}\right\}\,.\end{aligned}$$ This formula simultaneously captures the correct behaviour of $\Omega(p,k)$ in both soft limits $k \to 0$ and $k \to p$. It will be used in section \[sec:softAmpl\] below to compute the soft contributions to the reduced $2\to 2$ amplitude. Before doing that let us have a closer look at the $\epsilon$ expansion of the soft wavefunction we obtained.

We recall [@Caron-Huot:2017zfo] that all the negative powers of $\epsilon$ in (\[Well-1-ansatz-sym\]) cancel upon performing the sum over $n$, leading to a finite wavefunction at any loop order. While positive powers of $\epsilon$ in (\[Well-1-ansatz-sym\]) do play a role in the computation of the amplitude, the leading ${\mathcal{O}}({\epsilon}^0)$ have a special role: according to eq. (\[WhardTwod\]) it is precisely what must be subtracted from the full two-dimensional wavefunction to obtain the hard wavefunction ${{\Omega}_{\mathrm{h}}}^{({\rm 2d})}$. With this in mind, let us write down explicitly the leading terms in $\epsilon$ in the first few orders of the soft wavefunction in (\[Well-1-ansatz-sym\]):

\[Wsoft2dim\] $$\begin{aligned}
    {{{\Omega}^{({0})}}_{\mathrm{s}}}(p,k) \big|_{{\mathcal{O}}({\epsilon}^0)} &= 0, \\
    {{{\Omega}^{({1})}}_{\mathrm{s}}}(p,k) \big|_{{\mathcal{O}}({\epsilon}^0)} &= 
    \frac{{({C_A}-{{\mathbf{T}}_t^2})}}{2} \log {\left(\frac{k^2(p-k)^2}{(p^2)^2}\right)}, \\
    {{{\Omega}^{({2})}}_{\mathrm{s}}}(p,k) \big|_{{\mathcal{O}}({\epsilon}^0)} &= 
    \frac{{({C_A}-{{\mathbf{T}}_t^2})}^2}{4} \log^2 {\left(\frac{k^2(p-k)^2}{(p^2)^2}\right)}, \\
    {{{\Omega}^{({3})}}_{\mathrm{s}}}(p,k) \big|_{{\mathcal{O}}({\epsilon}^0)} &= 
    \frac{{({C_A}-{{\mathbf{T}}_t^2})}^3}{8} \log^3 {\left(\frac{k^2(p-k)^2}{(p^2)^2}\right)} 
    + \frac{{(2{C_A}-{{\mathbf{T}}_t^2})}{({C_A}-{{\mathbf{T}}_t^2})}^2}{2} \zeta_3, \\ \nn
    {{{\Omega}^{({4})}}_{\mathrm{s}}}(p,k) \big|_{{\mathcal{O}}({\epsilon}^0)} &= 
    \frac{{({C_A}-{{\mathbf{T}}_t^2})}^4}{16} \log^4 {\left(\frac{k^2(p-k)^2}{(p^2)^2}\right)} \\ 
    & +\, {(2{C_A}-{{\mathbf{T}}_t^2})}{({C_A}-{{\mathbf{T}}_t^2})}^3 \log {\left(\frac{k^2(p-k)^2}{(p^2)^2}\right)} \zeta_3, \\ \nn
    {{{\Omega}^{({5})}}_{\mathrm{s}}}(p,k) \big|_{{\mathcal{O}}({\epsilon}^0)} &= 
    \frac{{({C_A}-{{\mathbf{T}}_t^2})}^5}{32} \log^5 {\left(\frac{k^2(p-k)^2}{(p^2)^2}\right)} \\ 
    &+\, 
    \frac{{(2{C_A}-{{\mathbf{T}}_t^2})}{({C_A}-{{\mathbf{T}}_t^2})}^4}{4} \left[ 5 \log^2 {\left(\frac{k^2(p-k)^2}{(p^2)^2}\right)} \zeta_3 
    + 6 \zeta_5 \right]\,.\end{aligned}$$

In fact, these terms exponentiate and can be resummed into the following all-order expression using (\[OmegaEven\]) for $\epsilon=0$, yielding \[eq:wffullsoftresummed\] [\_]{}(p,k) |\_[(\^0)]{} = \^\^[ [([C\_A]{}-[\_t\^2]{})]{}]{}, with $x = L\, {\alpha_s}/\pi$.

Soft contributions to the $2\to2$ amplitude\[sec:softAmpl\]
-----------------------------------------------------------

Next, let us consider the soft contribution to the reduced $2\to 2$ scattering amplitude $\hat{\cal M}$. It is straighforward to insert [eq. ]{} into [eq. ]{}, perform the last integration and derive the $\ell$-th order contribution to the amplitude. In particular, given the symmetrised form of [eq. ]{}, the last integration can be done with the integration measure $[{\rm D}k]$ in [eq. ]{}, i.e. avoiding the need to introduce a cut-off as in ref. [@Caron-Huot:2017zfo]. After some arrangement we get \[MellReggeSoft-S\] \_[[NLL,s]{}]{}\^[(+,)]{} &=& i [([C\_A]{}-[\_t\^2]{})]{}\^[-1]{} \_[n=1]{}\^[-1]{} (-1)\^[n+1]{}\
&& \_[m=0]{}\^[n-2]{} [\_[s-u]{}\^2]{} [\^]{}, where the functions $B_n({\epsilon})$ and ${\hat B_{n}}({\epsilon})$ have been defined respectively in [eqs.  and ]{}, and we have introduced B\_[n]{}() = e\^[\_E]{} .

The coefficients ${\hat{\mathcal{M}}}_{{\rm NLL,s}}^{(+,\ell)}$ in (\[MellReggeSoft-S\]) are of course polynomial in the colour factors. For illustration, we expand [eq. ]{} to the first few orders in perturbation theory, obtaining

$$\begin{aligned}
    \label{eq:MsoftExpanded_1} 
    {\hat{\mathcal{M}}}_{\rm NLL,s}^{(1)} &= 
    i\pi {B_{0}}  \bigg\{\frac{1}{2{\epsilon}} \bigg\}  {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}, \\ 
    {\hat{\mathcal{M}}}_{\rm NLL,s}^{(2)} &= 
    i\pi \frac{{B_{0}}^2}{2} \bigg\{ \frac{{C_2}}{4 {\epsilon}^2} \bigg\} {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}, \\  
    {\hat{\mathcal{M}}}_{\rm NLL,s}^{(3)} &= 
    i\pi \frac{{B_{0}}^3}{3!} 
    \bigg\{ {C_2}^2 \left( \frac{1}{8 {\epsilon}^3} 
    - \frac{11\zeta_3}{4}\right) 
    - {C_1}{C_2}\frac{3\zeta_3}{4}\bigg\} {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}, \\  
    {\hat{\mathcal{M}}}_{\rm NLL,s}^{(4)} &= 
    i\pi \frac{{B_{0}}^4}{4!}
    \bigg\{ {C_2}^3 \left( \frac{1}{16 {\epsilon}^4} 
    + \frac{\zeta_3}{8{\epsilon}} + \frac{3\zeta_4}{16}\right) 
    + {C_1}{C_2}^2 \left( -\frac{\zeta_3}{8 {\epsilon}} 
    - \frac{3\zeta_4}{16}\right) \bigg\} {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}, \\  \nn
    {\hat{\mathcal{M}}}_{\rm NLL,s}^{(5)} &= 
    i\pi \frac{{B_{0}}^5}{5!}
    \bigg\{ {C_2}^4 \left( \frac{1}{32 {\epsilon}^5} 
    + \frac{\zeta_3}{16{\epsilon}^2} + \frac{3\zeta_4}{32 {\epsilon}} 
    - \frac{717\zeta_5}{16} \right) \\  
    &\hspace{1.0cm} +\, {C_1}{C_2}^3 \left( -\frac{\zeta_3}{16 {\epsilon}^2} 
    - \frac{3\zeta_4}{32 {\epsilon}} 
    - \frac{27\zeta_5}{16} \right) \bigg\} {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}, \\ \nn  
    {\hat{\mathcal{M}}}_{\rm NLL,s}^{(6)} &= 
    i\pi \frac{{B_{0}}^6}{6!}
    \bigg\{ {C_1}^2 {C_2}^3  \bigg( -\frac{39 \zeta_3^2}{16} \bigg) 
    + {C_1}{C_2}^4 \bigg(-\frac{\zeta_3}{32 {\epsilon}^3} 
    - \frac{3\zeta_4}{64 {\epsilon}^2}  
    - \frac{3\zeta_5}{32 {\epsilon}}
    - \frac{963\zeta_3^2}{32} 
    + \frac{5 \zeta_6}{32} \bigg) \\ 
    &\hspace{1.0cm}
    + {C_2}^5 \bigg(\frac{1}{64 {\epsilon}^6} 
    + \frac{\zeta_3}{32 {\epsilon}^3} 
    + \frac{3 \zeta_4}{64 {\epsilon}^2}  
    + \frac{3 \zeta_5}{32 {\epsilon}} 
    - \frac{2879 \zeta_3^2}{32} 
    + \frac{5 \zeta_6}{32} \bigg)
    \bigg\} {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}, \\  \nn
    {\hat{\mathcal{M}}}_{\rm NLL,s}^{(7)} &= 
    i\pi \frac{{B_{0}}^7}{7!}
    \bigg\{ {C_1}^2 {C_2}^4 \bigg( \frac{\zeta_3^2}{32 {\epsilon}} 
    + \frac{3\zeta_3 \zeta_4}{32} \bigg) 
    + {C_1}{C_2}^5 \bigg(
    - \frac{\zeta_3}{64 {\epsilon}^4} 
    - \frac{3\zeta_4}{128 {\epsilon}^3} 
    - \frac{3\zeta_5}{64 {\epsilon}^2}\\ \nn
    &\hspace{1.0cm}
    -\, \frac{3\zeta_3^2}{64 {\epsilon}} 
    - \frac{5 \zeta_6}{64 {\epsilon}} 
    - \frac{9 \zeta_3 \zeta_4}{64}
    - \frac{729 \zeta_7}{64} \bigg)  
    + {C_2}^6 \bigg( \frac{1}{128 {\epsilon}^7}
    + \frac{\zeta_3}{64 {\epsilon}^4} 
    + \frac{3 \zeta_4}{128 {\epsilon}^3} 
    + \frac{3 \zeta_5}{64 {\epsilon}^2} \\ 
    &\hspace{1.0cm}
    +\, \frac{\zeta_3^2}{64 {\epsilon}} 
    + \frac{5 \zeta_6}{64 {\epsilon}}  
    + \frac{3 \zeta_3 \zeta_4}{64} 
    - \frac{90711 \zeta_7}{64} \bigg)
    \bigg\} {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}, \\   \nn
    \label{eq:MsoftExpanded_8}
    {\hat{\mathcal{M}}}_{\rm NLL,s}^{(8)} &= 
    i\pi \frac{{B_{0}}^8}{8!}
    \bigg\{ {C_1}^2 {C_2}^5 \bigg( 
    \frac{\zeta_3^2}{64 {\epsilon}^2} 
    + \frac{3\zeta_3 \zeta_4}{64 {\epsilon}}
    - \frac{1341 \zeta_3 \zeta_5}{32} 
    + \frac{21\zeta_8}{512} \bigg)
    + {C_1}{C_2}^6 \bigg(
    -\frac{\zeta_3}{128 {\epsilon}^5} \\ \nn
    &\hspace{1.0cm}
    -\, \frac{3\zeta_4}{256 {\epsilon}^4}    
    - \frac{3\zeta_5}{128 {\epsilon}^3} 
    - \frac{3 \zeta_3^2}{128 {\epsilon}^2} 
    - \frac{5 \zeta_6}{128 {\epsilon}^2} 
    - \frac{9 \zeta_3 \zeta_4}{128 {\epsilon}} 
    - \frac{9\zeta_7}{128 {\epsilon}} 
    - \frac{96777 \zeta_3 \zeta_5}{64} \\ \nn
    &\hspace{1.0cm} 
    -\, \frac{189 \zeta_8}{1024} \bigg)
    + {C_2}^7 \bigg( \frac{1}{256 {\epsilon}^8} 
    + \frac{\zeta_3}{128 {\epsilon}^5} 
    + \frac{3\zeta_4}{256 {\epsilon}^4}
    + \frac{3\zeta_5}{128 {\epsilon}^3} 
    + \frac{\zeta_3^2}{128 {\epsilon}^2}
    + \frac{5 \zeta_6}{128 {\epsilon}^2}  \\ 
    &\hspace{1.0cm} 
    + \frac{3 \zeta_3 \zeta_4}{128 {\epsilon}} 
    + \frac{9 \zeta_7}{128 {\epsilon}} 
    - \frac{483837 \zeta_3 \zeta_5}{64} 
    + \frac{147 \zeta_8}{1024} \bigg)
    \bigg\} {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}\, ,
    \renewcommand{{C_1}}{(2{C_A}-{{\mathbf{T}}_t^2})}
    \renewcommand{{C_2}}{({C_A}-{{\mathbf{T}}_t^2})}\end{aligned}$$

where we used the shorthand notation for the colour factors, $C_1={C_1}$ and $C_2={C_2}$. We note that the expansion coefficients display uniform transcendental weight (where, as usual $1/\epsilon$ 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 $\zeta_2$ (or $\zeta_2$ times other zeta values, e.g. $\zeta_2\zeta_3$ at weight 5, etc.) factors do not appear in [eqs. –]{} ($\zeta_2$ terms would be present if we were to expand the factor factor $B_0^{\ell}({\epsilon})$). 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 coefficients ${\hat{\mathcal{M}}}_{{\rm NLL,s}}^{(+,\ell)}$ 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. [@Caron-Huot:2017zfo], so let us compare [eq. ]{} above with the result obtained in [@Caron-Huot:2017zfo]. There we proved that the singular terms of the reduced amplitude admit a simplified form \[MellReggeSoft-0\] \_[[NLL,s\_[simpl.]{}]{}]{}\^[(+,)]{} = [C\_2]{}\^[-1]{} ( 1- [B\_[-1]{}]{}() )\^[-1]{} [\_[s-u]{}\^2]{} [\^]{}. The latter, however, differs from the original soft amplitude obtained from [eq. ]{} starting at ${\cal O}({\epsilon}^{0})$ (compare eqs. (3.13) and (3.15) of [@Caron-Huot:2017zfo]). A nice feature of [eq. ]{} is that the loop functions in the amplitude at order $\ell$ do not depend on the index $\ell$, apart from the factor $\left({B_0}({\epsilon})\right)^{\ell}/\ell!$, and this allows one to easily resum [eq. ]{} to all orders in perturbation theory, obtaining an expression for the integrated soft amplitude ${\hat{\mathcal{M}}}_{{\rm NLL,s_{\rm simpl.}}}$: $$\begin{aligned}
 \label{MellReggeSoft-res}
\begin{split}
    {\hat{\mathcal{M}}}_{{\rm NLL,s_{\rm simpl.}}} =
    \frac{i\pi}{L {C_2}} \Bigg\{ &
    \Big( e^{\frac{B_0}{2{\epsilon}} {C_2}x}-1 \Big) \frac{B_{-1}({\epsilon})}{B_0({\epsilon})}  \\
    &\quad\times \left( 1- {\hat B_{-1}}({\epsilon}) \frac{{C_1}}{{C_2}} \right)^{-1}
    \bigg\} \, {{\mathbf{T}}_{s-u}^2}\, {{\mathcal{M}}^{\mathrm{(tree)}}}+ {\mathcal{O}}({\epsilon}^0),
\end{split}\end{aligned}$$ with $x =  L \,{\alpha_s}/\pi $ (see also eq. (3.18) of [@Caron-Huot:2017zfo]). 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 ${\cal O}(\epsilon^0)$ corrections to the amplitude defined in our current symmetric scheme, [eq. ]{}. To this end we consider the coefficients defined as the finite part of the *difference* between those in soft amplitude, eq. (\[MellReggeSoft-S\]), and in its simplified version, eq. (\[MellReggeSoft-0\]): $$\begin{aligned}
\label{deltadef}
\begin{split}
    {\hat{\mathcal{M}}}_{{\rm NLL,s}}^{(+,\ell)} - 
    {\hat{\mathcal{M}}}_{{\rm NLL,s_{\rm simpl.}}}^{(+,\ell)}
    &\equiv \, i \pi \,  
    \hat\Delta^{(+,\ell)}_{\rm NLL} \, {{\mathbf{T}}_{s-u}^2}\, {{\mathcal{M}}^{\mathrm{(tree)}}}+{\cal O}(\epsilon^1)\\ 
    &\equiv i \pi \, 
    \delta^{(\ell)} \, {C_2}^{\ell-1} \,{{\mathbf{T}}_{s-u}^2}\, {{\mathcal{M}}^{\mathrm{(tree)}}}+{\cal O}(\epsilon^1)\,.
\end{split}\end{aligned}$$ After some arrangement the coefficients $\hat\Delta^{(+,\ell)}_{\rm NLL}$ can be put into the form \[Delta\_s3\] \^[(+,)]{}\_[NLL]{} &=& [C\_2]{}\^[-1]{} { \_[n=0]{}\^[-1]{} (-1)\^[n+1]{}\
&& ( 1- [B\_[n-1]{}]{}() )\^[-1]{} \_[m=0]{}\^[n-2]{}}, where we discarded powers of $B_{0}(\epsilon)$, which do not affect the finite terms. From eq. (\[Delta\_s3\]) the coefficients $\delta^{(\ell)}$ of (\[deltadef\]) can be determined explicitly in terms of odd $\zeta$ numbers and the ratio of colour factors $r=\frac{{C_1}}{{C_2}}$. They are found to exponentiate in terms of the following rescaled odd $\zeta$ numbers: \[rescaled\_zeta\] \_[1+2 n]{} = \_[1+2n]{} (1+ ), such that the sum: $$\begin{aligned}
\label{summing_soft_finite}
\begin{split}
    \sum_{\ell=1}^{\infty} \frac{X^{\ell}}{\ell !}\delta^{(\ell)} 
    =\,&1-\exp\left(\sum_{n=1}^{\infty} X^{2n+1} \tilde\zeta_{2n+1}\right)
    \\=\,& 1-e^{-\gamma_E r X}
    \frac{\Gamma\Big(1 - X\Big)}{\Gamma\Big(1 + X\Big)} 
    \frac{\Big[\Gamma\Big(1+\frac{X}{2}\Big)\Big]^{2-r}}{
    \Big[\Gamma\Big(1 -  \frac{X}{2}\Big)\big]^{2-r}} 
\end{split}\end{aligned}$$ with $X\equiv {C_2}x$ and $x = L \,{\alpha_s}/\pi $, where we used 2 \_[n=1]{}\^ =-2 x \_E +((1-x))-((x+1)) . We conclude that the series $\hat\Delta^{(+,\ell)}_{\rm NLL}$ exponentiate to $$\begin{gathered}
\label{DeltaResum2} 
    \hat\Delta^{(+)}_{\rm NLL} =
    \frac{1}{L {C_2}} \Bigg[ 1- e^{-\gamma_E {C_1}\, x } 
    \frac{\Gamma\Big(1 - {C_2}x\Big)}{\Gamma\Big(1 + {C_2}x\Big)} \\ 
    \times \,
    \left(\frac{\Gamma\Big(1+{C_2}\frac{x}{2}\Big)}{\Gamma
    \Big(1 - {C_2}\frac{x}{2}\Big)}\right)^{-\frac{{\mathbf T}_t^2}{C_A-{\mathbf T}_t^2}}
    \Bigg].\end{gathered}$$ Using now the fact that that the simplified amplitude ${\hat{\mathcal{M}}}_{{\rm NLL,s_{\rm simpl.}}}^{(+,\ell)} $ in [eq. ]{} exponentiates independently, see [eq. ]{}, we obtain $$\begin{gathered}
\label{ReducedAmpNLLresum2} 
    {\hat{\mathcal{M}}}_{\rm NLL,s} = \\
    i\pi \,\Bigg\{ \frac{e^{\frac{B_0}{2{\epsilon}} {C_2}x}-1}{L {C_2}}
    \frac{B_{-1}({\epsilon})}{B_0({\epsilon})} 
    \left( 1- {\hat B_{-1}}({\epsilon}) \frac{{C_1}}{{C_2}} \right)^{-1}\,+ 
    \hat\Delta^{(+)}_{\rm NLL} \Bigg\} {{\mathbf{T}}_{s-u}^2}\, {{\mathcal{M}}^{\mathrm{(tree)}}}\\ 
    =
    i\pi \,\Bigg\{ \frac{e^{\frac{B_0}{2{\epsilon}} {C_2}x}-1}{L {C_2}}
    \left( 1- \frac{C_A}{{C_2}} R({\epsilon})\right)^{-1}\,+ 
    \hat\Delta^{(+)}_{\rm NLL} \Bigg\} {{\mathbf{T}}_{s-u}^2}\, {{\mathcal{M}}^{\mathrm{(tree)}}}\, ,\end{gathered}$$ where in the second line we expressed the amplitude in terms of the function $R({\epsilon}) ={{B_0}({\epsilon})}/{{B_{-1}}({\epsilon})} -1 $ of [eq. ]{}. Writing the reduced amplitude as in the second line of [eq. ]{} makes it easier to extract the infrared-renormalized amplitude from the reduced amplitude, as we will see in section \[hardFdef\]. Writing explicitly the factor $\hat\Delta^{(+)}_{\rm NLL} $, the reduced amplitude reads $$\begin{gathered}
\label{ReducedAmpNLLresum2B} 
    {\hat{\mathcal{M}}}_{\rm NLL,s} 
    = \frac{i\pi}{L {C_2}} \Bigg\{\left( e^{\frac{B_0}{2{\epsilon}} {C_2}x}-1 \right)  
    \left( 1- \frac{C_A}{{C_2}} R({\epsilon})\right)^{-1}\,+ 1 \\
    - e^{-\gamma_E {C_1}\, x} \,\,
    \frac{\Gamma\Big(1 - {C_2}x\Big)}{\Gamma\Big(1 + {C_2}x\Big)}
    \left(\frac{\Gamma\Big(1+{C_2}\frac{x}{2}\Big)}{\Gamma
    \Big(1 - {C_2}\frac{x}{2}\Big)}\right)^{-\frac{{\mathbf T}_t^2}{C_A-{\mathbf T}_t^2}}
    \Bigg\} {{\mathbf{T}}_{s-u}^2}\, {{\mathcal{M}}^{\mathrm{(tree)}}}.\end{gathered}$$ Of course, upon expansion (\[ReducedAmpNLLresum2B\]) yields back the coefficients of (\[MellReggeSoft-S\]) we listed in (\[eq:MsoftExpanded\_1\]) through (\[eq:MsoftExpanded\_8\]). Having at hand a resummed expression we can gain further insight on number-theoretical features of the expansion coefficients in [eqs. –]{}. We already know based on the derivation above that the $\hat\Delta^{(+)}_{\rm NLL}$ component in (\[ReducedAmpNLLresum2\]) gives rise to odd zeta values only. It then transpires that the sole origin of even ones is the function $R(\epsilon)$ in the first term. Further number-theoretical features will be discussed in section \[amplitude\], 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 ${\cal O}(\epsilon^0)$ terms is highly nontrivial, and it is an additional advantage of $k \leftrightarrow (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 {#hardFdef}
-----------------------------------------------------------------

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. ]{}. Following [eqs.  and ]{} we split the infrared-renormalized amplitude into a soft and a hard component: \[Hsplit\] \^[(+)]{}\_[NLL]{} = \^[(+)]{}\_[NLL,s]{} +\^[(+)]{}\_[NLL,h]{}. Then, from [eq. ]{} it follows that

\[getH2both\] $$\begin{aligned}
 \label{getH2-s} 
    {\mathcal{H}}^{(+)}_{\rm NLL,s} &=
    - e^{-\frac{x}{2{\epsilon}} {\mathbf{T}}_t^2} {\bf Z}^{(-)}_{\rm NLL} 
    \, e^{\frac{{B_{0}}({\epsilon})-1}{2{\epsilon}} \, x \, C_A} {\mathcal{M}}^{\rm (tree)}
    +e^{\frac{B_0({\epsilon})-1}{2{\epsilon}} \, x \, {\mathbf{T}}_t^2} {\cal \hat M}^{(+)}_{\rm NLL,s}, \\
    \label{getH2-h}
    {\mathcal{H}}^{(+)}_{\rm NLL,h} &= {\cal \hat M}^{(+)}_{\rm NLL,h}\,,\end{aligned}$$

where in (\[getH2-h\]) we neglected positive powers of $\epsilon$ originating in the expansion of $(1-B_0(\epsilon))$, using the fact that ${\cal \hat M}^{(+)}_{\rm NLL,h}$ is itself finite. Of course such a simplification cannot be applied to (\[getH2-s\]) where there is an interplay between positive powers of $\epsilon$ and negative ones. In section \[sec:softAmpl\] 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 $\epsilon\to 0$, according to [eq. ]{}. Inserting eqs. (\[Znll\]) and (\[ReducedAmpNLLresum2\]) into [eq. ]{} we get $$\begin{gathered}
\label{getH3}
    {\mathcal{H}}^{(+)}_{\rm NLL,s} = i\pi \Bigg\{ 
    - e^{\frac{{B_{0}}({\epsilon})-1}{2{\epsilon}} \, x \, C_A} 
    \bigg[\frac{ e^{\frac{x}{2{\epsilon}} ({C_{A}}-{{\mathbf{T}}_t^2})}}{L({C_{A}}-{{\mathbf{T}}_t^2})}
    \left( 1 - \frac{C_A}{C_A -{{\mathbf{T}}_t^2}} R({\epsilon}) \right)^{-1} \bigg]_{\rm poles} \\ 
    +e^{\frac{B_0({\epsilon})-1}{2{\epsilon}} \, x \, {\mathbf{T}}_t^2} 
    \Bigg[ \frac{e^{\frac{B_0}{2{\epsilon}} x {C_2}}-1}{L {C_2}}
    \left( 1- \frac{C_A}{{C_2}} R({\epsilon}) \right)^{-1}\,+ 
    \hat\Delta^{(+)}_{\rm NLL} \bigg]\Bigg\} {{\mathbf{T}}_{s-u}^2}\, {{\mathcal{M}}^{\mathrm{(tree)}}}\,,\end{gathered}$$ where we recall that $x = L \, {\alpha_s}/\pi$, and in the first line, corresponding to ${\bf Z}^{(-)}_{\rm NLL}$, we have dropped the $-1$ term in the numerator inside the square brackets, which does not generate any poles (see discussion following [eq. ]{}). This expression can be rearranged as follows: first of all, by collecting a factor $e^{\frac{{B_{0}}({\epsilon})-1}{2{\epsilon}} 
\, x \, C_A}$ we get $$\begin{gathered}
\label{getH4first}
    {\mathcal{H}}^{(+)}_{\rm NLL,s} = i\pi 
    e^{\frac{{B_{0}}({\epsilon})-1}{2{\epsilon}} \, x \, C_A}
    \Bigg\{-\bigg[\frac{ e^{\frac{x}{2{\epsilon}} ({C_{A}}-{{\mathbf{T}}_t^2})}}{L({C_{A}}-{{\mathbf{T}}_t^2})}
    \left( 1 - \frac{C_A}{C_A -{{\mathbf{T}}_t^2}} R({\epsilon}) \right)^{-1} \bigg]_{\rm poles} \\ 
    + \Bigg[ \frac{e^{\frac{x}{2{\epsilon}} {C_2}}
    -e^{\frac{1-B_0({\epsilon})}{2{\epsilon}} \, x \, {C_2}}}{L {C_2}}
    \left( 1- \frac{C_A}{{C_2}} R({\epsilon}) \right)^{-1} \\ 
    + e^{\frac{1-B_0({\epsilon})}{2{\epsilon}} \, x \, {C_2}}
    \hat\Delta^{(+)}_{\rm NLL} \bigg]\Bigg\} {{\mathbf{T}}_{s-u}^2}\, {{\mathcal{M}}^{\mathrm{(tree)}}}\,.\end{gathered}$$ We see at this point that the second line nicely cancel the poles from the first line. Furthermore, given that $1-B_0({\epsilon}) ={\cal O}({\epsilon}^2)$, see eq. (\[B0\]), and both $\big[ 1 - C_A/(C_A -{{\mathbf{T}}_t^2}) \, R({\epsilon}) \big]^{-1}
= 1 + {\mathcal{O}}({\epsilon}^3)$ and $\hat\Delta^{(+)}_{\rm NLL} = {\mathcal{O}}({\epsilon}^0)$, it is safe to set to one all exponentials containing the factor $1-B_0({\epsilon})$. We thus obtain $$\begin{gathered}
\label{getH4}
    {\mathcal{H}}^{(+)}_{\rm NLL,s} = i\pi 
    \Bigg\{ \bigg[\frac{ e^{\frac{x}{2{\epsilon}} ({C_{A}}-{{\mathbf{T}}_t^2})}-1}{L({C_{A}}-{{\mathbf{T}}_t^2})}
    \left( 1 - \frac{C_A}{C_A -{{\mathbf{T}}_t^2}} R({\epsilon}) \right)^{-1} \bigg]_{{\epsilon}^0} 
    + \hat\Delta^{(+)}_{\rm NLL} \Bigg\} {{\mathbf{T}}_{s-u}^2}\, {{\mathcal{M}}^{\mathrm{(tree)}}}\,,\end{gathered}$$ with $\hat\Delta^{(+)}_{\rm NLL}$ given in [eq. ]{}. For later reference we also expand [eq. ]{} to the first few orders in perturbation theory, obtaining (recall $C_2=C_A-{{\mathbf{T}}_t^2}$):

$$\begin{aligned}
    \label{eq:HsoftExpanded_1} 
    {\cal H}_{\rm NLL,s}^{(1)} &= 0, \\ 
    {\cal H}_{\rm NLL,s}^{(2)} &= 0, \\  
    {\cal H}_{\rm NLL,s}^{(3)} &= \frac{i\pi}{3!}  
    \bigg\{ - C_A {C_2}\frac{3\zeta_3}{4}
    -{C_2}^2 \frac{7\zeta_3}{2} \bigg\} {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}, \\  
    {\cal H}_{\rm NLL,s}^{(4)} &= \frac{i\pi}{4!}
    \bigg\{ - C_A {C_2}^2 \frac{3\zeta_4}{16} \bigg\} {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}, \\ 
    {\cal H}_{\rm NLL,s}^{(5)} &= \frac{i\pi}{5!}
    \bigg\{- C_A {C_2}^3 \frac{27\zeta_5}{16}
    -{C_2}^4 \frac{93\zeta_5}{2} \bigg\} {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}, \\ 
    {\cal H}_{\rm NLL,s}^{(6)} &= \frac{i\pi}{6!}
    \bigg\{ - C_A^2 {C_2}^3  \frac{39 \zeta_3^2}{16} 
    -C_A {C_2}^4 \bigg(
    \frac{1119\zeta_3^2}{32} 
    + \frac{5\zeta_6}{32} \bigg) 
    - {C_2}^5 \frac{245 \zeta_3^2}{2} \bigg\} {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}, \\  
    {\cal H}_{\rm NLL,s}^{(7)} &= \frac{i\pi}{7!}
    \bigg\{ C_A^2 {C_2}^4 \frac{3\zeta_3 \zeta_4}{32} 
    + C_A {C_2}^5 \bigg( 
    \frac{3\zeta_3 \zeta_4}{64}
    - \frac{729 \zeta_7}{64} \bigg)  
    - {C_2}^6 \frac{5715 \zeta_7}{4} \bigg\} {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}, \\   \nn
    \label{eq:HsoftExpanded_8}
    {\cal H}_{\rm NLL,s}^{(8)} &= \frac{i\pi}{8!}
    \bigg\{ C_A^2 {C_2}^5 \bigg(
    - \frac{1341 \zeta_3 \zeta_5}{32} 
    + \frac{21\zeta_8}{512} \bigg)
    + C_A {C_2}^6 \bigg(
    - \frac{102141 \zeta_3 \zeta_5}{64}  
    - \frac{105\zeta_8}{1024} \bigg) \\
    &\hspace{5.0cm} \,- {C_2}^7 \, 9114 \zeta_3 \zeta_5\bigg\} {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}\, .
\renewcommand{{C_1}}{(2{C_A}-{{\mathbf{T}}_t^2})}
\renewcommand{{C_2}}{({C_A}-{{\mathbf{T}}_t^2})}\end{aligned}$$

It is interesting note that $\zeta_n$ values with even $n$ originate solely from the expansion of the factor $R({\epsilon})$ in [eq. ]{}, while the expansion of the factor $\hat\Delta^{(+)}_{\rm NLL}$ generates only $\zeta_n$ values with odd $n$. The latter property of $\hat\Delta^{(+)}_{\rm NLL}$ makes this function compatible with the class of zeta values we will encounter considering the two-dimensional amplitude in section \[amplitude\].

In summary, according to (\[Hsplit\]) the infrared-renormalized amplitude is given as a sum of two terms: $\mathcal{H}_{\rm s}$, computed in this section using the soft approximation, plus $\mathcal{H}_{\rm h}$, which is identical to the hard part of the reduced amplitude (see [eq. ]{}). 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 {#2d-bfkl}
===========================================

As discussed in the introduction and in section \[chap:bfkl\], much of the complication of solving the 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 (\[OmegaSplitDef\]) subject to the condition (\[h\_definition\]): the vanishing of ${{\Omega}_{\mathrm{h}}}^{({\rm 2d})}(p,k)$ in the soft limits guaranties that the corresponding amplitude $\hat {\cal M}^{(+)}_{\rm NLL,h}\left(\frac{s}{-t}\right)$ of eq. (\[Mhard\]) is finite.

Our task in this section is therefore to compute ${{\Omega}_{\mathrm{h}}}^{({\rm 2d})}(p,k)$. We do so by iteratively applying the Hamiltonian of eqs. (\[Hdef1\]) and (\[Hamil\]), according to eq. (\[Hdef0\]). We keep the kinematics general, but in contrast to section \[chap:bfkl\], 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 k = k\_x + i k\_y, k’ = k\_x’ + i k\_y’ p = p\_x + i p\_y, where the real and imaginary parts are the components of the corresponding momenta and introduce new variables $z,w \in \mathbb{C}$ according to \[eq:zwdef\] = = . Since the wavefunction is a function of Lorentz scalars (i.e. squares of momenta) it will be symmetric under the exchange $z \leftrightarrow {\bar{z}}$ with ${\bar{z}}$ the complex conjugate of $z$. In particular, ${{\Omega}_{\mathrm{h}}}^{({\rm 2d})}(p,k)$ depends on the two ratios $$\label{k2zzb}
    \frac{k^2}{p^2}= \frac{z{\bar{z}}}{(1-z)(1-{\bar{z}})}\,,\quad\qquad 
    \frac{(p-k)^2}{p^2}= \frac{1}{(1-z)(1-{\bar{z}})}\,.$$ These relations also clarify that the symmetry under interchanging the two Reggeons, eq. (\[left-right-symmetry\]), corresponds to $z\to1/z$, and specifically, the two soft limits where one or the other Reggeon is soft correspond respectively to $z\to 0$ and $z\to \infty$. The limit $z\to 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 kernel reads \[eq:f2d\] p\^2 f(p,k,k’) (1-w)\^2 (1-[|[w]{}]{})\^2 K(w,[|[w]{}]{},z,[|[z]{}]{}), where \[eq:K2d\] K(w,[|[w]{}]{},z,[|[z]{}]{}) = = + + . Furthermore, in the limit ${\epsilon}\to 0$, $J(p,k)$ of eq.  becomes \[eq:j2d\] J(p,k) j(z,[|[z]{}]{}) 12 , and the measure reads \[eq:dw2d\] . Here, ${\mathrm{d}}^2 w \equiv {\mathrm{d}}{{\rm Re}}(w)\,{\mathrm{d}}{{\rm Im}}(w)$ where the real and imaginary part of $w$ are to be integrated from $-\infty$ to $+\infty$, in accordance with eq. .

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 ${{\Omega}_{\mathrm{h}}}^{({\rm 2d})}(p,k)= {{\Omega}_{\mathrm{2d}}}(z,{\bar{z}})$. We expand it as in eq. (\[OmegaEven\]), where we take $B_0(0)= 1$, i.e. \[OmegaEven\_hard\] [\_]{}(z,[|[z]{}]{}) &=& \_[=0]{}\^ ( )\^[+1]{} [[\^[()]{}]{}\_]{}(z,[|[z]{}]{}), where the coefficients of increasing orders are related by the action of the Hamiltonian according to eq. (\[Hdef0\]), which now reads: \[eq:Htdaction\] [[\^[()]{}]{}\_]{}(z,[|[z]{}]{}) = [\_]{}[[\^[([-1]{})]{}]{}\_]{}(z,[|[z]{}]{}), where [\_]{}(z,[|[z]{}]{}) = [C\_1]{}[\_]{}(z,[|[z]{}]{}) + [C\_2]{}[\_]{}(z,[|[z]{}]{}). Plugging in the above expressions we find the two parts of the Hamiltonian to be

$$\begin{aligned}
    \label{eq:hi2d} {{\hat{H}}_{\mathrm{2d,i}}}\Psi(z,{\bar{z}}) &= \frac{1}{4\pi} 
    \int {\mathrm{d}}^2 w K(w,{\bar{w}},z,{\bar{z}}) \left[ \Psi(w,{\bar{w}}) - \Psi(z,{\bar{z}}) \right], \\
    \label{eq:hm2d} {{\hat{H}}_{\mathrm{2d,m}}}\Psi(z,{\bar{z}}) &= j(z,{\bar{z}}) \Psi(z,{\bar{z}}),\end{aligned}$$

where ${{{\Omega}^{({0})}}_{\mathrm{2d}}}(z,{\bar{z}}) = {{\Omega}^{({0})}}(p,k) = 1$.

In the next section we proceed to solve for the wavefunction ${{\Omega}_{\mathrm{2d}}}$ by iterating the two-dimensional Hamiltonian .

The two-dimensional wavefunction {#sec:wf2d}
--------------------------------

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* (), which were first described by Brown in ref. [@Brown:2004ugm]. Since then, several applications of in computing scattering amplitudes have been found, in particular in the context of the high-energy limit, e.g. [@Pennington:2012zj; @Dixon:2012yy; @DelDuca:2013lma; @Dixon:2014voa; @DelDuca:2016lad; @DelDuca:2018hrv], and in the context of infared singularities in general kinematics [@Almelid:2017qju; @Dixon:2019lnw]. Here we will show that these functions also form a suitable basis for expressing the two-dimensional wavefunction ${{{\Omega}^{({\ell})}}_{\mathrm{2d}}}(z,{\bar{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 () of $z$ with of ${\bar{z}}$. We shall denote by ${\mathcal{L}}_\sigma(z,{\bar{z}})$ where $\sigma$ is a sequence of *letters*, typically zeros and ones.[^3] The *letters* are said to form an *alphabet*, $\{0,1\}$, and $\sigma$ is, by analogy, referred to as a *word*. The length of a word is often called the (transcendental) *weight* of the .

are the natural choice for the two-dimensional evolution, since $j(z,{\bar{z}})$ of eq.  belongs to this class, j(z,[|[z]{}]{}) = 12 \_[0]{}(z,[|[z]{}]{}) + \_[1]{}(z,[|[z]{}]{}), and the action of the Hamiltonian ${{\hat{H}}_{\mathrm{2d,i}}}$ preserves single-valuedness when acting on a single-valued function. This can be expected on general grounds: any complex pair $z,{\bar{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, single-valuedness 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 \[app:svhpls\].

As noted upon introducing the variables $z$ and $\bar{z}$ in (\[eq:zwdef\]), the two-dimensional wavefunction is symmetric under $z \leftrightarrow {\bar{z}}$. In addition, as mentioned following (\[k2zzb\]), owing to the symmetry upon interchanging the two Reggeons in [eq. ]{}, the wavefunction is invariant under simultaneously swapping $z \leftrightarrow 1/z$ and ${\bar{z}}\leftrightarrow 1/{\bar{z}}$. Both these symmetries are easily verified by looking at eqs.  and , where, for the latter symmetry, one changes the integration variables $w \to 1/w$, ${\bar{w}}\to 1/{\bar{w}}$. We will use these properties to simplify the iteration of the wavefunction as well as its results in section \[sec:asalphabet\].

The evolution of the wavefunction in strictly two transverse dimensions according to (\[eq:Htdaction\]) has the following basic characteristics. Firstly, iterating ${{\hat{H}}_{\mathrm{2d,m}}}$ amounts to multiplying by $j(z,{\bar{z}})$ and therefore evaluating shuffle products of . Secondly, each application of ${{\hat{H}}_{\mathrm{2d,i}}}$ adds one layer of integration such that ${{{\Omega}^{({\ell-1})}}_{\mathrm{2d}}}$ can be written as a linear combination of of weight $\ell-1$. A method to calculate the convolution in eq. (\[eq:hi2d\]) in terms of residues was described in chapter 6 of Ref. [@DelDuca:2018hrv]. 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 .

Suppose we wish to compute the action of a linear operator $\hat O$, which may involve integration, on a function $\Psi(z,{\bar{z}})$. Assume now that we find a differential operator $\Delta$, which is linear in logarithmic derivatives with respect to $z$ and ${\bar{z}}$, with the following properties:

\[eq:deltacrit\] $$\begin{aligned}
    \label{eq:deltacriti} i.\ &\Delta \text{ commutes with } \hat{O} \\
    \label{eq:deltacritii} ii.\ &\Delta \Psi \text{ is a pure function with a 
    weight that is lower than } \Psi~\text{by one unit}.\end{aligned}$$

Then, \[eq:deltadiffeq\] = , and we can compute $\hat{O} \left[ \Psi(z,{\bar{z}}) \right]$ by integrating the differential equation , 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 $\Delta$ is assumed to be linear in derivatives with respect to $z$ and ${\bar{z}}$, solving the differential equation amounts to computing a one-dimensional integral. This may be contrasted with the original integral in (\[eq:hi2d\]) which is two-dimensional. Given , solving this differential equation is straightforward, and the result remains in the class of , (see eq. ). The same applies for the class of : to solve the differential equation within this class, we first integrate its holomorphic part according to eq. , and subsequently recover the full result, depending on both $z$ and ${\bar{z}}$, by applying the single-valued map $\mathbf{s}$ defined in eq. . Having outlined the general approach let us see how it is implemented in practice to solve for the wavefunction in (\[eq:hi2d\]).

Let us start by considering the $\hat{O}$ in eq. (\[eq:deltadiffeq\]) to coincide with the two-dimensional Hamiltonian ${{\hat{H}}_{\mathrm{2d,i}}}$ (we will see below that the final procedure involves considering parts of the Hamiltonian ${{\hat{H}}_{\mathrm{2d,i}}}$ in turn). The most natural candidate for the operator $\Delta$ in eq. (\[eq:deltadiffeq\]) is $\Delta_1=z {{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}$, since condition (\[eq:deltacriti\]) is satisfied, as we now show.

For generic values of $w$ and $z$ one finds using eq.  \[eq:zddzKsym\] zK(w,[|[w]{}]{},z,[|[z]{}]{}) = -w K(w,[|[w]{}]{},z,[|[z]{}]{}) . This implies that $z {{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}$ commutes with the Hamiltonian, $$\begin{aligned}
    z{{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}\left[ {{\hat{H}}_{\mathrm{2d,i}}}\Psi(z,{\bar{z}}) \right] &= \frac{1}{4\pi} 
    \int {\mathrm{d}}^2 w \left\{ \left( -{{\frac{{\mathrm{d}}}{{\mathrm{d}}{w}}}}w K(w,{\bar{w}},z,{\bar{z}}) \right) \left[ \Psi(w,{\bar{w}}) - \Psi(z,{\bar{z}}) \right] \right. \nn \\
    &\hspace{45mm} \left. - K(w,{\bar{w}},z,{\bar{z}}) \left( z{{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}\Psi(z,{\bar{z}}) \right) \right\} \nn \\
    &= \frac{1}{4\pi} \int {\mathrm{d}}^2 w K(w,{\bar{w}},z,{\bar{z}}) \left[ w{{\frac{{\mathrm{d}}}{{\mathrm{d}}{w}}}}\Psi(w,{\bar{w}}) - z{{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}\Psi(z,{\bar{z}}) \right] \nn \\
    \label{eq:zddzhisym} &= {{\hat{H}}_{\mathrm{2d,i}}}\left[ z{{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}\Psi(z,{\bar{z}}) \right] \qquad \text{(for generic $w,z$)}.\end{aligned}$$ fulfilling condition (\[eq:deltacriti\]). However, some extra caution is needed here: the complex-conjugate pairs $w,{\bar{w}}$ and $z,{\bar{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 two-dimensional Poisson equation \[eq:2dpoisson\] \_w \_[[|[w]{}]{}]{} (w [|[w]{}]{}) = \^2(w) namely, by contributions of the form \[eq:cterms\] = \^2(w-c) with $c$ a complex number. The two-dimensional $\delta$ function in the above equations fixes both the real and the imaginary part of its argument such that \^2 w \^2(w-c) f(w,[|[w]{}]{}) = f(c,|[c]{}) for some function $f$, *cf.* the remark below eq. .

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. . Eq.  therefore correctly reads $$\begin{aligned}
    z{{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}K(w,{\bar{w}},z,{\bar{z}}) &\!= \!\left[ z{{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}K(w,{\bar{w}},z,{\bar{z}}) \right]_{\mathrm{reg}}\!\!
    + \left[ z{{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}K(w,{\bar{w}},z,{\bar{z}}) \right]_{\mathrm{con}}\nn \\
    &\!=\! -\left[ {{\frac{{\mathrm{d}}}{{\mathrm{d}}{w}}}}w K(w,{\bar{w}},z,{\bar{z}}) \right]_{\mathrm{reg}}\!\!
    + \left[ z{{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}K(w,{\bar{w}},z,{\bar{z}}) \right]_{\mathrm{con}}\nn \\
    &\!=\! -{{\frac{{\mathrm{d}}}{{\mathrm{d}}{w}}}}\left[ w K(w,{\bar{w}},z,{\bar{z}}) \right] 
    + \left[ {{\frac{{\mathrm{d}}}{{\mathrm{d}}{w}}}}w K(w,{\bar{w}},z,{\bar{z}}) \right]_{\mathrm{con}}\!\!
    + \left[ z{{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}K(w,{\bar{w}},z,{\bar{z}}) \right]_{\mathrm{con}}\end{aligned}$$ which modifies eq.  to give $$\begin{gathered}
    \label{eq:zddzhisym2} z{{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}\left[ {{\hat{H}}_{\mathrm{2d,i}}}\Psi(z,{\bar{z}}) \right] 
    = {{\hat{H}}_{\mathrm{2d,i}}}\left[ z{{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}\Psi(z,{\bar{z}}) \right] \\
    + \frac{1}{4\pi} \int {\mathrm{d}}^2 w \left\{ \left[ {{\frac{{\mathrm{d}}}{{\mathrm{d}}{w}}}}w K(w,{\bar{w}},z,{\bar{z}}) \right]_{\mathrm{con}}+ \left[ z{{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}K(w,{\bar{w}},z,{\bar{z}}) \right]_{\mathrm{con}}\right\} \\
    \times \left[ \Psi(w,{\bar{w}}) - \Psi(z,{\bar{z}}) \right].\end{gathered}$$ We shall continue to refer to the behaviour in eq.  as the commutativity of $z{{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}$ and ${{\hat{H}}_{\mathrm{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) $\delta$-function which makes the integral on the r.h.s. of eq.  easy to evaluate.

We will derive the explicit form of the contact terms towards the end of this section, at which point eq. (\[eq:zddzhisym2\]) will become directly usable for determining the action of ${{\hat{H}}_{\mathrm{2d,i}}}$ on the wavefunction $\Psi$. Before doing that, however, we turn our attention to condition (\[eq:deltacritii\]). Concretely in eq. (\[eq:zddzhisym2\]) the requirement is that $z\frac{d}{dz} \Psi$ should be a pure function of weight ones less than $\Psi$ itself. We find that the operator $z\frac{d}{dz}$, upon acting on any of the form ${\mathcal{L}}_{0,\sigma}(z,{\bar{z}})$, does indeed yield such a pure function, so eq. (\[eq:zddzhisym2\]) becomes: z = [\_]{} + \[eq:zddzL0\], where we have used eq. . On the other hand, $z{{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}$ does not have the same effect when acting on an ${\mathcal{L}}_{1,\sigma}(z,{\bar{z}})$, where one obtains instead z = [\_]{}+ , which does not fulfil the condition . One may be tempted to use $(1-z){{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}$ instead but, unfortunately, this operator does not commute with ${{\hat{H}}_{\mathrm{2d,i}}}$.

The solution is to first split the Hamiltonian ${{\hat{H}}_{\mathrm{2d,i}}}= {{\hat{H}}_{\mathrm{2d,i}_1}}+ {{\hat{H}}_{\mathrm{2d,i}_2}}$ with [\_[\_n]{}]{}(z,[|[z]{}]{}) = \^2 w K\_n(w,[|[w]{}]{},z,[|[z]{}]{}) and

$$\begin{aligned}
    K_1(w,{\bar{w}},z,{\bar{z}}) &= 
    \left(\frac{1}{w-z}-\frac{1}{w}\right) \frac{1}{{\bar{w}}-{\bar{z}}} \\
    K_2(w,{\bar{w}},z,{\bar{z}}) &= 
    \frac{1}{w-z} \left(\frac{1}{{\bar{w}}-{\bar{z}}}-\frac{1}{{\bar{w}}}\right)\end{aligned}$$

where $K_1(w,{\bar{w}},z,{\bar{z}}) + K_2(w,{\bar{w}},z,{\bar{z}}) 
= K(w,{\bar{w}},z,{\bar{z}})$, *cf.* eq. . This split is useful because it opens the possibility of identifying different differential operators $\Delta_i$ that commute with the separate components of the Hamiltonian ${{\hat{H}}_{\mathrm{2d,i}_1}}$ and ${{\hat{H}}_{\mathrm{2d,i}_2}}$, *and* yield a pure function when acting directly on ${\mathcal{L}}_{0,\sigma}(z,{\bar{z}})$ or on ${\mathcal{L}}_{1,\sigma}(z,{\bar{z}})$, thus simultaneously fulfilling both conditions in .

Regarding the commutation relations, condition , it is straightforward to verify that the following four relations hold, up to contact terms:

\[eq:Hicommut\] $$\begin{aligned}
    \left[ z{{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}, {{\hat{H}}_{\mathrm{2d,i}_1}}\right] &= {\text{(contact terms)}}, & 
    \left[ z(1-z){{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}, {{\hat{H}}_{\mathrm{2d,i}_1}}\right] &= {\text{(contact terms)}}\label{eq:Hi1ops}, \\
    \left[ z{{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}, {{\hat{H}}_{\mathrm{2d,i}_2}}\right] &= {\text{(contact terms)}}, & 
    \left[ (1-z){{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}, {{\hat{H}}_{\mathrm{2d,i}_2}}\right] &= {\text{(contact terms)}}. \label{eq:Hi2ops} \end{aligned}$$

Let us therefore define the following three differential operators: \[Deltai\] \_i = f\_i(z) f\_i(z) = {

[ll]{} z& i=1\
1-z& i=2\
z(1-z) & i=3

., and show that we can arrange the wavefunction, which is a linear combination of ${\mathcal{L}}_{0,\sigma}(z,{\bar{z}})$ and ${\mathcal{L}}_{1,\sigma}(z,{\bar{z}})$, such that condition (\[eq:deltacritii\]) would also be fulfilled.

To this end, let us first note that upon acting on ${\mathcal{L}}_{0,\sigma}(z,{\bar{z}})$ with either of the two parts of the Hamiltonian we have (using ): $$\label{eq:hi1diffeq2} z{{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}\left[ {{\hat{H}}_{\mathrm{2d,i}_n}}{\mathcal{L}}_{0,\sigma}(z,{\bar{z}})  \right] 
    = {{\hat{H}}_{\mathrm{2d,i}_n}}\left[ {\mathcal{L}}_\sigma(z,{\bar{z}}) \right] + {\text{(contact terms)}}\,,$$ just as in (\[eq:zddzL0\]). Thus, the remaining challenge is to handle terms containing ${\mathcal{L}}_{1,\sigma}(z,{\bar{z}})$; this is where the additional flexibility of splitting the Hamiltonian pays off. Let us consider first the simplest case of ${{\hat{H}}_{\mathrm{2d,i}_2}}$ where we obtain \[eq:hi2diffeq\] (1-z) = [\_[\_2]{}]{}+ . Now ${{\hat{H}}_{\mathrm{2d,i}_2}}\Psi$ can be readily integrated for any $\Psi$ using (\[eq:hi1diffeq2\]) and (\[eq:hi2diffeq\]). Turning to consider ${{\hat{H}}_{\mathrm{2d,i}_1}}$, let us write \_[1,]{}(z,[|[z]{}]{}) = (\_[1,]{}(z,[|[z]{}]{}) + \_[0,]{}(z,[|[z]{}]{})) - \_[0,]{}(z,[|[z]{}]{}) \[eq:L1split\] and use the linearity of the Hamiltonian to act with it on $({\mathcal{L}}_{1,\sigma} + {\mathcal{L}}_{0,\sigma})$ and $(-{\mathcal{L}}_{0,\sigma})$ separately. We may now apply respectively the differential operators $\Delta_3$ and $\Delta_1$ of to these terms. With eq.  and one can easily verify that they produce the desired pure functions of lower weight in accordance with : $$\begin{aligned}
    \label{eq:hi1diffeq1} z(1-z){{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}\left[ {{\hat{H}}_{\mathrm{2d,i}_1}}\left( {\mathcal{L}}_{0,\sigma}(z,{\bar{z}}) 
    + {\mathcal{L}}_{1,\sigma}(z,{\bar{z}}) \right) \right] 
    &= {{\hat{H}}_{\mathrm{2d,i}_1}}\left[ {\mathcal{L}}_\sigma(z,{\bar{z}}) \right] + {\text{(contact terms)}}\end{aligned}$$ Using (\[eq:hi1diffeq1\]) along with (\[eq:hi1diffeq2\]) we see that also ${{\hat{H}}_{\mathrm{2d,i}_1}}\Psi$ can be integrated for any $\Psi$. 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 (\[eq:deltacrit\]).

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.  for general $\Delta_i = f_i(z) {{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}$ and the two parts of the split Hamiltonian, $$\begin{gathered}
    \label{eq:fddzhinsym} f_i(z) {{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}\left[ {{\hat{H}}_{\mathrm{2d,i}_n}}\Psi(z,{\bar{z}}) \right] 
    = {{\hat{H}}_{\mathrm{2d,i}_n}}\left[ f_i(z) {{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}\Psi(z,{\bar{z}}) \right] \\
    + \frac{1}{4\pi} \int {\mathrm{d}}^2 w \left\{ \left[ {{\frac{{\mathrm{d}}}{{\mathrm{d}}{w}}}}f_i(w) K_n(w,{\bar{w}},z,{\bar{z}}) \right]_{\mathrm{con}}+ \left[ f_i(z) {{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}K_n(w,{\bar{w}},z,{\bar{z}}) \right]_{\mathrm{con}}\right\} \\
    \times \left[ \Psi(w,{\bar{w}}) - \Psi(z,{\bar{z}}) \right]\end{gathered}$$ where, according to eqs. (\[eq:hi1diffeq2\]), (\[eq:hi2diffeq\]) and (\[eq:hi1diffeq1\]), the relevant combinations of $i$ and $n$ are

$$\begin{aligned}
    n = 1 \quad &\longrightarrow \quad i = 1 \text{ or } 3 \\
    n = 2 \quad &\longrightarrow \quad i = 1 \text{ or } 2\,.\end{aligned}$$

In computing the contact terms in (\[eq:fddzhinsym\]) we note that the $f_i(z)$  are functions of $z$ only whilst being independent of the complex conjugate ${\bar{z}}$. According to eq.  this implies that \[eq:pulloutf\] \_ = f\_i(w) \_. for $n=1,2$, and thus (\[eq:fddzhinsym\]) becomes: $$\begin{gathered}
    \label{eq:fddzhinsym_simpl} f_i(z) {{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}\left[ {{\hat{H}}_{\mathrm{2d,i}_n}}\Psi(z,{\bar{z}}) \right] 
    = {{\hat{H}}_{\mathrm{2d,i}_n}}\left[ f_i(z) {{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}\Psi(z,{\bar{z}}) \right] \\
    + \frac{1}{4\pi} \int {\mathrm{d}}^2 w \left\{ \left[f_i(w) {{\frac{{\mathrm{d}}}{{\mathrm{d}}{w}}}}K_n(w,{\bar{w}},z,{\bar{z}}) \right]_{\mathrm{con}}+ \left[ f_i(z) {{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}K_n(w,{\bar{w}},z,{\bar{z}}) \right]_{\mathrm{con}}\right\} \\
    \times \left[ \Psi(w,{\bar{w}}) - \Psi(z,{\bar{z}}) \right]\,.\end{gathered}$$ Consequently, we only have to consider the following four derivatives,

$$\begin{aligned}
    \label{eq:ddwk1ct} \left[ {{\frac{{\mathrm{d}}}{{\mathrm{d}}{w}}}}K_1(w,{\bar{w}},z,{\bar{z}}) \right]_{\mathrm{con}}&= \pi \left[ \delta^2(w-z) - \delta^2(w-\infty) \right] \frac{z}{w(w-z)} \\
    \label{eq:ddzk1ct} \left[ {{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}K_1(w,{\bar{w}},z,{\bar{z}}) \right]_{\mathrm{con}}&= -\pi \delta^2(z-w) \frac{z}{w(w-z)} \\[3mm]
    \label{eq:ddwk2ct} \left[ {{\frac{{\mathrm{d}}}{{\mathrm{d}}{w}}}}K_2(w,{\bar{w}},z,{\bar{z}}) \right]_{\mathrm{con}}&= \pi \left[ \delta^2(w-z) - \delta^2(w) \right] \frac{1}{w-z} \\  
    \label{eq:ddzk2ct} \left[ {{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}K_2(w,{\bar{w}},z,{\bar{z}}) \right]_{\mathrm{con}}&= -\pi \delta^2(z-w) \frac{1}{w-z},\end{aligned}$$

where in eqs.  and we have dropped terms proportional to $\delta^2(z)$, restricting our calculation to $z \neq 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.  the terms proportional to $\delta^2(w-z) = \delta^2(z-w)$ in eqs. – cancel identically, so the remaining contact-term contributions are only at $w=\infty$ for $K_1$ and at $w=0$ for $K_2$. Using the corresponding $\delta$ functions to turn the integrals over $w$ in (\[eq:fddzhinsym\_simpl\] ) into evaluation of limits at infinity and at zero respectively we finally obtain: $$\begin{aligned}
    \label{eq:fddzhinsym_H1} f_i(z) {{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}\left[ {{\hat{H}}_{\mathrm{2d,i}_1}}\Psi(z,{\bar{z}}) \right] 
    &= {{\hat{H}}_{\mathrm{2d,i}_1}}\left[ f_i(z) {{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}\Psi(z,{\bar{z}}) \right] 
    -\frac{1}{4} \lim_{w\to\infty}  \frac{z f_i(w)}{w(w-z)}  
    \left[ \Psi(w,{\bar{w}}) - \Psi(z,{\bar{z}}) \right]\,, \\
    \label{eq:fddzhinsym_H2} f_i(z) {{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}\left[ {{\hat{H}}_{\mathrm{2d,i}_2}}\Psi(z,{\bar{z}}) \right] 
    &= {{\hat{H}}_{\mathrm{2d,i}_2}}\left[ f_i(z) {{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}\Psi(z,{\bar{z}}) \right] 
    - \frac{1}{4}  \lim_{w\to 0}  \frac{f_i(w) }{w-z}  
    \left[ \Psi(w,{\bar{w}}) - \Psi(z,{\bar{z}}) \right]\,.\end{aligned}$$ This equations will be used in the next section to determine the wavefunction.

Differential equations and an iterative solution for the wavefunction {#sec:diffeq}
---------------------------------------------------------------------

Finding the differential equations is now simply a matter of compiling together the results of the previous section. Starting with the easiest case, $\Delta_1 {{\hat{H}}_{\mathrm{2d,i}_n}}{\mathcal{L}}_{0,\sigma}$, we notice that with $f_1(w) = w$ both the $w\to \infty$ limit in eq. (\[eq:fddzhinsym\_H1\]) and the $w\to 0$ limit in eq. (\[eq:fddzhinsym\_H2\]) vanish, and thus there are no contributions from contact terms in either of these cases. Dividing by $f_1(z) = z$ to arrive at [\_[\_n]{}]{}\_[0,]{}(z,[|[z]{}]{}) = . \[eq:diffeq0\] Next consider the case $\Delta_2 {{\hat{H}}_{\mathrm{2d,i}_2}}{\mathcal{L}}_{1,\sigma}$, corresponding to eq. . Here $f_2(w) = 1-w$ and eq. (\[eq:fddzhinsym\_H2\]) yields \[eq:diffeq1b\] [\_[\_2]{}]{}\_[1,]{}(z,[|[z]{}]{}) = - 14 . where we have divided by $f_2(z) = 1-z$ and used the shorthand $[\ldots]_{w,{\bar{w}}\rightarrow 0}$ to denote the $w,{\bar{w}}\rightarrow 0$ limit of the functions inside the square brackets. This term can, in fact, be dropped as it always contains a single whose indices feature (at least) one “1” and, thus, is equal to zero in the limit.

The last case, $\Delta_i {{\hat{H}}_{\mathrm{2d,i}_1}}{\mathcal{L}}_{1,\sigma}$, is governed by eqs.  and , using the wavefunction split of eq. (\[eq:L1split\]). Considering in turn the action of eq. (\[eq:fddzhinsym\_H1\]) on $({\mathcal{L}}_{1,\sigma}(z,{\bar{z}}) +{\mathcal{L}}_{0,\sigma}(z,{\bar{z}}))$ with $f_3(w) 
= w(1-w)$ and on $(- {\mathcal{L}}_{0,\sigma}(z,{\bar{z}}))$ with $f_1(w) = w$, we derive two separate equations, which we then combine using the linearity of operators ${{\hat{H}}_{\mathrm{2d,i}_1}}$ and ${{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}$ to obtain $$\begin{gathered}
 \label{eq:diffeq1a}
    {{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}{{\hat{H}}_{\mathrm{2d,i}_1}}{\mathcal{L}}_{1,\sigma}(z,{\bar{z}}) 
    = \frac{{{\hat{H}}_{\mathrm{2d,i}_1}}{\mathcal{L}}_{\sigma}(z,{\bar{z}})}{1-z} \\
    - \frac14 \frac{{\mathcal{L}}_{0,\sigma}(z,{\bar{z}}) 
    + {\mathcal{L}}_{1,\sigma}(z,{\bar{z}}) - [{\mathcal{L}}_{0,\sigma}(w,{\bar{w}}) 
    + {\mathcal{L}}_{1,\sigma}(w,{\bar{w}})]_{w,{\bar{w}}\rightarrow \infty}}{1-z} \end{gathered}$$ with $[\ldots]_{w,{\bar{w}}\rightarrow \infty}$ the $w,{\bar{w}}\rightarrow \infty$ limit of the functions inside the square brackets. Taking this limit requires some careful analytic continuation of the relevant to ensure that $w$ and ${\bar{w}}$ stay complex-conjugate as they approach infinity.

Because the Hamiltonian ${{\hat{H}}_{\mathrm{2d,i}}}$ and its components ${{\hat{H}}_{\mathrm{2d,i}_n}}$ are linear operators one can sum up the above equations – and recombine ${{\hat{H}}_{\mathrm{2d,i}_1}}+ {{\hat{H}}_{\mathrm{2d,i}_2}}\to {{\hat{H}}_{\mathrm{2d,i}}}$ obtaining more compact expressions:

\[diffeq\] $$\begin{aligned}
    {{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}{{\hat{H}}_{\mathrm{2d,i}}}{\mathcal{L}}_{0,\sigma}(z,{\bar{z}}) 
    &= \frac{{{\hat{H}}_{\mathrm{2d,i}}}{\mathcal{L}}_{\sigma}(z,{\bar{z}})}{z}\,, 
    \label{eq:diffeq0new} \\
    {{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}{{\hat{H}}_{\mathrm{2d,i}}}{\mathcal{L}}_{1,\sigma}(z,{\bar{z}}) 
    &= \frac{{{\hat{H}}_{\mathrm{2d,i}}}{\mathcal{L}}_{\sigma}(z,{\bar{z}})}{1-z} 
    - \frac14 \frac{{\mathcal{L}}_{1,\sigma}(z,{\bar{z}})}{z} \nn \\
    &\hspace{-5mm} - \frac14 \frac{{\mathcal{L}}_{0,\sigma}(z,{\bar{z}}) 
    + 2{\mathcal{L}}_{1,\sigma}(z,{\bar{z}}) - [{\mathcal{L}}_{0,\sigma}(w,{\bar{w}}) 
    + {\mathcal{L}}_{1,\sigma}(w,{\bar{w}})]_{w,{\bar{w}}\rightarrow \infty}}{1-z}\,. 
    \label{eq:diffeq1new}\end{aligned}$$

These differential equations compactly represent the action of the Hamiltonian $\hat{H}_{{\rm 2d},{\rm i}}$ according to eq. (\[eq:hi2d\]). 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$ and ${\bar{z}}$ — a small detour is necessary to recover the action of ${{\hat{H}}_{\mathrm{2d,i}}}$ on : we take the holomorphic part of a given , integrate it w.r.t. $z$ according to the differential equations in (\[diffeq\]), and then reconstruct the functional dependence on ${\bar{z}}$ by requiring the result be single-valued. This ultimately amounts to simply replacing \_0\^z t \_[0,]{}(z,[|[z]{}]{}) \_0\^z t \_[1,]{}(z,[|[z]{}]{}) . For more details on this procedure see appendix \[app:holomorphicpart\].

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=z{\bar{z}}$ tend to zero. For small $z,{\bar{z}}$, only with all-zero indices can give non-zero contributions; these correspond to powers of logarithms: $${\mathcal{L}}_{\vec{0}_n}(z,{\bar{z}}) = \frac{\log^n(z{\bar{z}})}{n!} 
    \qquad \text{with} \qquad \vec{0}_n = 
    \underbrace{0,\dots,0}_{n\ \text{zeros}} \,.$$ In eq.  we calculated the action of the small-$k$ (or soft) Hamiltonian ${{\hat{H}}_{\mathrm{s}}}$ on powers of $\xi = (k^2/p^2)^{-{\epsilon}}$. The action of ${{\hat{H}}_{\mathrm{i}}}$ in the soft limit can be isolated by looking at the coefficient of $2{C_A}- {{\mathbf{T}}_t^2}$ and thus is $${{\hat{H}}_{\mathrm{i}}}|_{\text{soft}} \left( \frac{k^2}{p^2} \right)^{-m{\epsilon}} 
    = \frac{\hat B_m({\epsilon})}{2{\epsilon}} 
    \left( \frac{k^2}{p^2} \right)^{-(m+1){\epsilon}} \label{eq:Hik2p2}$$ where $\hat B_m({\epsilon})$ is given in eq. (\[bubblehat\]). Expanding both sides in ${\epsilon}$ and matching powers of $\delta=m\epsilon$ in the limit ${\epsilon}\rightarrow 0$ lets us extract the action of ${{\hat{H}}_{\mathrm{i}}}$ in the soft limit on any given power of $\log (k^2/p^2) = \log (z{\bar{z}})$. For reference, we find

$$\begin{aligned}
    {{\hat{H}}_{\mathrm{i}}}|_{\text{soft}} {\mathcal{L}}_0(z,{\bar{z}}) 
    &= \mathcal{O}({\epsilon}) \label{HiL0} \\
    {{\hat{H}}_{\mathrm{i}}}|_{\text{soft}} {\mathcal{L}}_{0,0}(z,{\bar{z}}) 
    &= \zeta_3 + \mathcal{O}({\epsilon}) \\
    {{\hat{H}}_{\mathrm{i}}}|_{\text{soft}} {\mathcal{L}}_{0,0,0}(z,{\bar{z}}) 
    &= \zeta_3 {\mathcal{L}}_0(z,{\bar{z}}) + \mathcal{O}({\epsilon}) \\
    {{\hat{H}}_{\mathrm{i}}}|_{\text{soft}} {\mathcal{L}}_{0,0,0,0}(z,{\bar{z}}) 
    &= \zeta_3 {\mathcal{L}}_{0,0}(z,{\bar{z}}) + \zeta_5 + \mathcal{O}({\epsilon}) \\
    {{\hat{H}}_{\mathrm{i}}}|_{\text{soft}} {\mathcal{L}}_{0,0,0,0,0}(z,{\bar{z}}) 
    &= \zeta_3 {\mathcal{L}}_{0,0,0}(z,{\bar{z}}) + \zeta_5 {\mathcal{L}}_0(z,{\bar{z}}) 
    + \mathcal{O}({\epsilon}) \label{HiL00000} \end{aligned}$$

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 expanding ${\hat{B}_{m}}({\epsilon})$ which is a product of gamma functions.

We can now calculate the action of ${{\hat{H}}_{\mathrm{2d,i}}}$ on any by iteratively solving the differential equations and , starting from the lowest-weight functions, ${\mathcal{L}}_0$ and ${\mathcal{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.  numerically and compare to the analytical result for a number of randomly generated pairs $z,{\bar{z}}$. Specifically, with $w = w_1 + i w_2$ and $z = z_1 + i z_2$ the action of ${{\hat{H}}_{\mathrm{2d,i}}}$ can be written $$\begin{gathered}
    {{\hat{H}}_{\mathrm{2d,i}}}\Psi(z,{\bar{z}}) = \frac{1}{2\pi} \int_{-\infty}^\infty 
    {\mathrm{d}}w_1 \int_{-\infty}^\infty {\mathrm{d}}w_2 \frac{w_1 z_1 
    + w_2 z_2}{(w_1^2 + w_2^2)((w_1 - z_1)^2 + (w_2 - z_2)^2)} \\
    \times \left[ \Psi(w_1 + i w_2,w_1 - i w_2) 
    - \Psi(z_1 + i z_2,z_1 - i z_2) \right].\end{gathered}$$ where $\Psi(z,{\bar{z}})$ is a (linear combination of) (s). This type of integral is readily evaluated numerically in e.g. `Mathematica`.

For the wavefunction up to weight four we find

$$\begin{aligned}
    {{{\Omega}^{({1})}}_{\mathrm{2d}}} &= \frac{1}{2} C_2 \left({\mathcal{L}}_0+2 {\mathcal{L}}_1\right) \label{eq:wtd1} \\
    {{{\Omega}^{({2})}}_{\mathrm{2d}}} &= \frac{1}{2} C_2^2 \left({\mathcal{L}}_{0,0}+2 {\mathcal{L}}_{0,1}
    +2 {\mathcal{L}}_{1,0}+4 {\mathcal{L}}_{1,1}\right)+\frac{1}{4} C_1 C_2 \left(-{\mathcal{L}}_{0,1}
    -{\mathcal{L}}_{1,0}-2 {\mathcal{L}}_{1,1}\right) \\
    {{{\Omega}^{({3})}}_{\mathrm{2d}}} &=
    \frac{3}{4} C_2^3 \left({\mathcal{L}}_{0,0,0}+2 {\mathcal{L}}_{0,0,1}
    +2 {\mathcal{L}}_{0,1,0}+4 {\mathcal{L}}_{0,1,1}+2 {\mathcal{L}}_{1,0,0}
    +4 {\mathcal{L}}_{1,0,1}+4 {\mathcal{L}}_{1,1,0}
    +8 {\mathcal{L}}_{1,1,1}\right)
    \nn \\  &\hspace{4mm}+
    \frac{1}{4} C_1 C_2^2 \left(2 \zeta_3 -2 {\mathcal{L}}_{0,0,1}
    -3 {\mathcal{L}}_{0,1,0}-7 {\mathcal{L}}_{0,1,1}-2 {\mathcal{L}}_{1,0,0}
    -7 {\mathcal{L}}_{1,0,1}-7 {\mathcal{L}}_{1,1,0}
    -14 {\mathcal{L}}_{1,1,1}\right)
    \nn \\&\hspace{4mm}
    +\frac{1}{16} C_1^2 C_2 \left({\mathcal{L}}_{0,0,1}
    +2 {\mathcal{L}}_{0,1,0}+4 {\mathcal{L}}_{0,1,1} 
    +{\mathcal{L}}_{1,0,0}+4 {\mathcal{L}}_{1,0,1}
    +4 {\mathcal{L}}_{1,1,0}+8 {\mathcal{L}}_{1,1,1}\right)\\
    {{{\Omega}^{({4})}}_{\mathrm{2d}}} &= 
    \frac{3}{2} C_2^4 \left({\mathcal{L}}_{0,0,0,0}
    +2 {\mathcal{L}}_{0,0,0,1}+2 {\mathcal{L}}_{0,0,1,0}+4 {\mathcal{L}}_{0,0,1,1}
    +2 {\mathcal{L}}_{0,1,0,0}+4 {\mathcal{L}}_{0,1,0,1}\right. \nn \\
    &\hspace{4mm} \left. +4 {\mathcal{L}}_{0,1,1,0}+8 {\mathcal{L}}_{0,1,1,1}
    +2 {\mathcal{L}}_{1,0,0,0}+4 {\mathcal{L}}_{1,0,0,1}+4 {\mathcal{L}}_{1,0,1,0}
    +8 {\mathcal{L}}_{1,0,1,1}\right. \nn \\
    &\hspace{4mm} \left. +4 {\mathcal{L}}_{1,1,0,0}+8 {\mathcal{L}}_{1,1,0,1}
    +8 {\mathcal{L}}_{1,1,1,0}+16 {\mathcal{L}}_{1,1,1,1}\right)
    \nn \\
    &\hspace{4mm}
    +\frac{1}{8} C_1 C_2^3 \left(-9 {\mathcal{L}}_{0,0,0,1}
    -14 {\mathcal{L}}_{0,0,1,0}-34 {\mathcal{L}}_{0,0,1,1}-14 {\mathcal{L}}_{0,1,0,0}
    -42 {\mathcal{L}}_{0,1,0,1}\right. \nn \\
    &\hspace{4mm} \left. -44 {\mathcal{L}}_{0,1,1,0}-92 {\mathcal{L}}_{0,1,1,1}
    -9 {\mathcal{L}}_{1,0,0,0}-34 {\mathcal{L}}_{1,0,0,1}-42 {\mathcal{L}}_{1,0,1,0}
    -92 {\mathcal{L}}_{1,0,1,1}\right. \nn \\
    &\hspace{4mm} \left. -34 {\mathcal{L}}_{1,1,0,0}-92 {\mathcal{L}}_{1,1,0,1}
    -92 {\mathcal{L}}_{1,1,1,0}-184 {\mathcal{L}}_{1,1,1,1}+8 {\mathcal{L}}_0 \zeta_3
    +28 {\mathcal{L}}_1 \zeta_3\right) \nn \\
    &\hspace{4mm} +\frac{1}{32} C_1^2 C_2^2 \left(7 {\mathcal{L}}_{0,0,0,1}
    +15 {\mathcal{L}}_{0,0,1,0}+34 {\mathcal{L}}_{0,0,1,1}+15 {\mathcal{L}}_{0,1,0,0}
    +56 {\mathcal{L}}_{0,1,0,1}\right. \nn \\
    &\hspace{4mm} \left. +56 {\mathcal{L}}_{0,1,1,0}+116 {\mathcal{L}}_{0,1,1,1}
    +7 {\mathcal{L}}_{1,0,0,0}+40 {\mathcal{L}}_{1,0,0,1}+56 {\mathcal{L}}_{1,0,1,0}\right. \nn \\
    &\hspace{4mm} \left. +116 {\mathcal{L}}_{1,0,1,1}+34 {\mathcal{L}}_{1,1,0,0}
    +116 {\mathcal{L}}_{1,1,0,1}+116 {\mathcal{L}}_{1,1,1,0}+232 {\mathcal{L}}_{1,1,1,1}
    -44 {\mathcal{L}}_1 \zeta_3\right) \nn \\
    &\hspace{4mm}
    +\frac{1}{64} C_1^3 C_2 \left(-{\mathcal{L}}_{0,0,0,1}
    -3 {\mathcal{L}}_{0,0,1,0}-6 {\mathcal{L}}_{0,0,1,1}-3 {\mathcal{L}}_{0,1,0,0}\right. \nn \\
    &\hspace{4mm} \left. -12 {\mathcal{L}}_{0,1,0,1}-12 {\mathcal{L}}_{0,1,1,0}
    -24 {\mathcal{L}}_{0,1,1,1}-{\mathcal{L}}_{1,0,0,0}-8 {\mathcal{L}}_{1,0,0,1}
    -12 {\mathcal{L}}_{1,0,1,0}\right. \nn \\
    &\hspace{4mm} \left. -24 {\mathcal{L}}_{1,0,1,1}-6 {\mathcal{L}}_{1,1,0,0}
    -24 {\mathcal{L}}_{1,1,0,1}-24 {\mathcal{L}}_{1,1,1,0}-48 {\mathcal{L}}_{1,1,1,1}
    +12 {\mathcal{L}}_1 \zeta_3\right) 
    \label{eq:wtd4}\end{aligned}$$

where we introduced the notation $C_1 = 2{C_A}- {{\mathbf{T}}_t^2}$, $C_2 = {C_A}- {{\mathbf{T}}_t^2}$ and wrote ${{{\Omega}^{({\ell})}}_{\mathrm{2d}}} \equiv {{{\Omega}^{({\ell})}}_{\mathrm{2d}}}(z,{\bar{z}})$ and ${\mathcal{L}}_\sigma \equiv {\mathcal{L}}_\sigma(z,{\bar{z}})$ for brevity. Further results up to weight 14 can be found in the ancillary file `2Reggeon-wavefunction-L01-Basis.txt`.

Interestingly, a new type of transcendental number appears for the first time in the twelve-loop wavefunction — a so-called multiple zeta value (). While it is no surprise that 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 number-theoretical implications.

are the values of evaluated at special points, typically their branch points $z=1$ or $z \to \infty$, for example[^4] $H_{0,0,0,0,1,0,0,1}(1) = 
H_{5,3}(1) = \zeta_{5,3}$. It turns out that only cover a subset of all when evaluated at $z = {\bar{z}}= 1$ or $z,{\bar{z}}\to \infty$ and we refer to this subset as single-valued multiple zeta values. They are discussed in detail in refs. [@Schnetz:2013hqa; @Brown:2013gia] where the authors show that, up to weight ten, the algebra of single-valued is generated by ordinary (odd) zeta numbers $\zeta_n$. At weight eleven, however, a new type of number appears, alongside the expected $\zeta_{11}$. We shall call it[^5] $g_{5,3,3}$ and it is defined by g\_[5,3,3]{} = -47 \_2\^3 \_5 +65 \_2\^2 \_7 + 45 \_2 \_9 + \_[5,3,3]{}, \[eq:g533\] where $\zeta_{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,{\bar{w}}\to \infty$ limit in eq. . The former are generated by expanding gamma functions, *cf.* eq.  with eq. , and can therefore contribute only *single* (ordinary) zeta numbers. The value of the large-$w,{\bar{w}}$ limit instead does generally involve (single-valued) *multiple* zeta values. We note that it is guaranteed to multiply the weight-one ${\mathcal{L}}_1(z,{\bar{z}})$ which is generated by the denominator, $1-z$, upon integrating the differential equation . Being the sole source of (single-valued) , we conclude that such zeta values of weight $w$ can only occur starting at the next loop order, i.e. $\ell=w+1$. Specifically, this explains why $g_{5,3,3}$, which is weight 11, cannot appear at loop orders $\ell < 12$. Indeed, we find that $g_{5,3,3}$ is accompanied by ${\mathcal{L}}_1$ in the twelve-loop wavefunction: $$\begin{gathered}
    {{{\Omega}^{({12})}}_{\mathrm{2d}}}(z,{\bar{z}}) \supset \frac{1}{80} \left( 
    \frac{88653 C_2^2 C_1^{10}}{2048}
    -\frac{1021171 C_2^3 C_1^9}{4096}
    -\frac{3517129 C_2^4 C_1^8}{1024} \right. \\
    +\frac{43378313 C_2^5 C_1^7}{1024}
    -\frac{5951395 C_2^6 C_1^6}{32}
    +\frac{1583033 C_2^7 C_1^5}{4} \\
    \left. -\frac{6320709 C_2^8 C_1^4}{16}
    + 135513 C_2^9 C_1^3 \right) \times g_{5,3,3} \, {\mathcal{L}}_1(z,{\bar{z}})\,.\end{gathered}$$ According to ref. [@Brown:2013gia] (*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 \[eq:g553\] g\_[5,5,3]{} = 10 \_2\^2 \_9 + \_2 \_[11]{} + 5 \_5 \_[5,3]{} + \_[5,5,3]{} and \[eq:g733\] g\_[7,3,3]{} = - \_2\^3 \_7 + \_2\^2 \_9 + \_2 \_[11]{} + 6 \_5 \_[5,3]{} + \_[7,3,3]{} we observe that the 14-loop wavefunction contains the term $$\begin{gathered}
    {{{\Omega}^{({14})}}_{\mathrm{2d}}}(z,{\bar{z}}) \supset \frac{1}{2240}\left( 
    - \frac{132291047 C_2^2 C_1^{12}}{20480} 
    + \frac{7701138629 C_2^3 C_1^{11}}{183500800} \right. \\
    - \frac{21177619993 C_2^4 C_1^{10}}{81920} 
    - \frac{141869475599 C_2^5 C_1^9}{40960} 
    + \frac{144180124197 C_2^6 C_1^8}{4096} \\
    - \frac{1550199662073 C_2^7 C_1^7}{10240} 
    + \frac{941115705999 C_2^8 C_1^6}{2560} 
    - \frac{41630406511 C_2^9 C_1^5}{80} \\ \left. 
    + \frac{15828500247 C_2^{10} C_1^4}{40} 
    - 120229353 C_2^{11} C_1^3 \right) \times g_{5,5,3} \, {\mathcal{L}}_1(z,{\bar{z}})\end{gathered}$$ as well as $$\begin{gathered}
    {{{\Omega}^{({14})}}_{\mathrm{2d}}}(z,{\bar{z}}) \supset \frac{1}{896} \left( 
    \frac{557319 C_2^2 C_1^{12}}{256} 
    - \frac{296956417 C_2^3 C_1^{11}}{16384} \right. \\
    - \frac{3811324785 C_2^4 C_1^{10}}{16384} 
    + \frac{36358896425 C_2^5 C_1^9}{8192} 
    - \frac{125984665967 C_2^6 C_1^8}{4096} \\
    + \frac{241764230539 C_2^7 C_1^7}{2048} 
    - \frac{139303244409 C_2^8 C_1^6}{512} 
    + \frac{11897473261 C_2^9 C_1^5}{32} \\ \left. 
    - \frac{2180551359 C_2^{10} C_1^4}{8} 
    + 79134813 C_2^{11} C_1^3 \right) \times g_{7,3,3} \, {\mathcal{L}}_1(z,{\bar{z}}).\end{gathered}$$ The observed term $g_{5,3,3}\, {\mathcal{L}}_1(z,{\bar{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 . The non-zero coefficients of $g_{5,5,3}\, {\mathcal{L}}_1(z,{\bar{z}})$ and $g_{7,3,3}\, {\mathcal{L}}_1(z,{\bar{z}})$ at 14 loops may be seen as hint that indeed *all* single-valued 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 when interpreting our results for the amplitude in section \[sec:finiteamp\].

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 {#sec:asalphabet}
------------------------

Throughout this paper we have tried to exploit the symmetries of the 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 \[sec:wf2d\], in two dimensions, the wavefunction is invariant under two transformations: complex conjugation and inversion of the arguments. The latter, i.e. the fact that ${{\Omega}_{\mathrm{2d}}}(z,{\bar{z}}) = {{\Omega}_{\mathrm{2d}}}(1/z,1/{\bar{z}})$, corresponds to eq. (\[left-right-symmetry\]), i.e. to swapping the two reggeons, and was used, for example, to identify the two soft limits in section \[soft\]. In the present context, it inspired us to introduce a new alphabet for , as we now explain. Instead of 0 and 1, corresponding to integration over ${\mathrm{d}}\log z$ and ${\mathrm{d}}\log (1-z)$, respectively, we shall use $a$ and $s$. They are associated with integration over ${\mathrm{d}}\log z$ and ${\mathrm{d}}\log z/(1-z)^2$ and thus behave antisymmetrically and symmetrically, respectively, under $z \to 1/z$. In particular $$\label{ElDef}
    {\mathcal{L}}_s(z,{\bar{z}})=\log \frac{z{\bar{z}}}{(1-z)^2(1-{\bar{z}})^2}
    \qquad\Longrightarrow\qquad 
    {\mathcal{L}}_s(1/z,1/{\bar{z}})={\mathcal{L}}_s(z,{\bar{z}})\,.$$ The leading-order wavefunction simplifies to ${{{\Omega}^{({1})}}_{\mathrm{2d}}} = \frac{1}{2} C_2 {\mathcal{L}}_s(z,{\bar{z}})$, and at higher orders, the $z \to 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,{\bar{z}}) = {\mathcal{L}}_s(z,{\bar{z}})/2$ of eq.  and hence the action of ${{\hat{H}}_{\mathrm{2d,m}}}$ in eq. , which now amounts to shuffling an $s$ into the indices of the function it acts on (and multiplying by a $\frac12$), for example [\_]{}\_[a,s,a,s]{}(z,[|[z]{}]{}) = 12 \_[s,a,s,a,s]{}(z,[|[z]{}]{}) + \_[a,s,s,a,s]{}(z,[|[z]{}]{}) + \_[a,s,a,s,s]{}(z,[|[z]{}]{}). The action of ${{\hat{H}}_{\mathrm{2d,i}}}$ has a much richer and more complicated structure. However, we notice that at symbol level, i.e. keeping only the highest-weight , it simply amounts to replacing $s \to ss - aa$ and multiplying by $-\frac14$, for example, $$\begin{gathered}
    \label{eq:hisymbollevel} {{\hat{H}}_{\mathrm{2d,i}}}{\mathcal{L}}_{a,s,a,s}(z,{\bar{z}}) 
    = -\frac14 \Big( {\mathcal{L}}_{a,s,s,a,s}(z,{\bar{z}}) - {\mathcal{L}}_{a,a,a,a,s}(z,{\bar{z}}) 
    + {\mathcal{L}}_{a,s,a,s,s}(z,{\bar{z}}) - {\mathcal{L}}_{a,s,a,a,a}(z,{\bar{z}}) \Big) \\
    + \Sigma_{\text{sub}}\,,\end{gathered}$$ where $\Sigma_{\text{sub}}$ contains products of subleading-weight and zeta numbers, i.e. terms like ${\mathcal{L}}_\sigma \zeta_{n_1} \cdots \zeta_{n_m}$ with $|\sigma| + n_1 + \dots + n_m = 5$ and $|\sigma| < 5$ in the above example. This replacement rule can be derived from the differential equations and , as we now explain.

To this end, let us consider the two cases ${{\hat{H}}_{\mathrm{2d,i}}}{\mathcal{L}}_{a,\sigma}$ and ${{\hat{H}}_{\mathrm{2d,i}}}{\mathcal{L}}_{s,\sigma}$ in turn. Considering the former, due to the equivalence of the letters 0 and $a$, eq.  immediately gives the action on ${\mathcal{L}}_{a,\sigma}$ \[eq:diffeqaleading\] [\_]{}\_[a,]{}(z,[|[z]{}]{}) = . The simple recursive nature of this equation implies that ${{\hat{H}}_{\mathrm{2d,i}}}$ does not affect the $a$ indices of a and can, at most, generate subleading terms $\Sigma_{\text{sub}}$ through integration constants, *cf.* eqs. –.

The action on ${\mathcal{L}}_{s,\sigma}$ can be written as $$\begin{aligned}
    {{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}{{\hat{H}}_{\mathrm{2d,i}}}{\mathcal{L}}_{s,\sigma}(z,{\bar{z}}) 
    &= {{\frac{{\mathrm{d}}}{{\mathrm{d}}{z}}}}{{\hat{H}}_{\mathrm{2d,i}}}\left[ {\mathcal{L}}_{0,\sigma}(z,{\bar{z}}) 
    + 2{\mathcal{L}}_{1,\sigma}(z,{\bar{z}}) \right] \nn\\
    &= \frac{1+z}{z(1-z)} {{\hat{H}}_{\mathrm{2d,i}}}{\mathcal{L}}_\sigma(z,{\bar{z}}) \nn \\
    &\hspace{10mm} - \frac12 \left( \frac{{\mathcal{L}}_{0,\sigma}(z,{\bar{z}}) 
        + 2{\mathcal{L}}_{1,\sigma}(z,{\bar{z}})}{1-z} 
    + \frac{{\mathcal{L}}_{1,\sigma}(z,{\bar{z}})}{z} \right) 
    + \Sigma_{\text{sub}} \nn \\
    &= \frac{1+z}{z(1-z)} {{\hat{H}}_{\mathrm{2d,i}}}{\mathcal{L}}_\sigma(z,{\bar{z}}) \nn \\
    &\hspace{10mm} - \frac14 \frac{1+z}{z(1-z)} 
    \left( {\mathcal{L}}_{0,\sigma}(z,{\bar{z}}) + 2{\mathcal{L}}_{1,\sigma}(z,{\bar{z}}) \right) 
    + \frac14 \frac{{\mathcal{L}}_{0,\sigma}(z,{\bar{z}})}{z} + \Sigma_{\text{sub}} \nn \\
    &= \frac{1+z}{z(1-z)} {{\hat{H}}_{\mathrm{2d,i}}}{\mathcal{L}}_\sigma(z,{\bar{z}}) \nn \\
    &\hspace{10mm} - \frac14 \left( \frac{1+z}{z(1-z)} {\mathcal{L}}_{s,\sigma}(z,{\bar{z}}) 
    - \frac{{\mathcal{L}}_{a,\sigma}(z,{\bar{z}})}{z} \right) + \Sigma_{\text{sub}}, \label{eq:diffeqsleading}\end{aligned}$$ where at each step we have used $\Sigma_{\text{sub}}$ to collect subleading terms into. The first term in the final expression is again an inert term, like the one encountered in eq. . The following term however, creates two leading-weight terms which, upon integration, yield $-\frac14 ({\mathcal{L}}_{s,s,\sigma} 
- {\mathcal{L}}_{a,a,\sigma})$ and hence confirm the pattern described above eq. . Note that by the recursive nature of the differential equation this applies (separately) to *every* letter $s$ in the word $(s,\sigma)$, not just the first one (see e.g. eq. (\[eq:hisymbollevel\])).

In the following we show that it is possible to unravel the recursive definition of ${{\hat{H}}_{\mathrm{2d,i}}}$ beyond symbol level. The $\Sigma_{\text{sub}}$ terms in the above equations are generated by two independent and additive sources: the $w,{\bar{w}}\to \infty$ limit in eq.  and the constants of integration as shown in eqs. –. Let us denote them $\Sigma_{\text{sub}(\infty)}$ and $\Sigma_{\text{sub}(0)}$, respectively, with their sum equalling $\Sigma_{\text{sub}}$. Empirically we observe that $\Sigma_{\text{sub}(0)}$ follows a simple pattern when using the $\{a,s\}$ alphabet: [\_]{}\_[w\_1,…,w\_[-1]{}]{}(z,[|[z]{}]{}) = \_ + \_[()]{} + \_[j 3 ]{}\^\_j \_[w\_1,…,w\_[-j]{}]{}(z,[|[z]{}]{}). \[eq:subtermszero\] with $\Sigma_{\text{lead}}$ now the leading-weight governed by eq. . $\Sigma_{\text{sub}(\infty)}$ in turn can be summarised by $$\begin{gathered}
    {{\hat{H}}_{\mathrm{i}}}{\mathcal{L}}_{w_1,\dots,w_{\ell-1}}(z,{\bar{z}}) 
    = \Sigma_{\text{lead}} + \Sigma_{\text{sub}(0)} 
    + \frac18 \sum_{j=1}^{\ell-1}  \left( {\mathcal{L}}_{w_1,\dots,w_j}(z,{\bar{z}}) 
    -{\mathcal{L}}_{w_1,\dots,w_{j-1},a}(z,{\bar{z}}) \right) \\
    \times \left[ {\mathcal{L}}_{a,w_{j+1},\dots,w_{\ell-1}}(z,{\bar{z}}) 
    + {\mathcal{L}}_{s,w_{j+1},\dots,w_{\ell-1}}(z,{\bar{z}}) \right]_{z,{\bar{z}}\to \infty}. 
    \label{eq:subtermsinf}\end{gathered}$$ In both these equations the final term in the sum needs to be interpreted with care: in eq. (\[eq:subtermszero\]), for $j=\ell$ one obtains ${\mathcal{L}}_{w_1,\ldots,w_0}\equiv1$ and in eq. (\[eq:subtermsinf\]) for $j=\ell-1$ one obtains in the second factor ${\mathcal{L}}_{a,w_{\ell},\dots,w_{\ell-1}}(z,{\bar{z}}) 
+ {\mathcal{L}}_{s,w_{\ell},\dots,w_{\ell-1}}(z,{\bar{z}}) \equiv
{\mathcal{L}}_{a}(z,{\bar{z}})  + {\mathcal{L}}_{s}(z,{\bar{z}})$. Observe that in eq. (\[eq:subtermsinf\]) is a necessary yet not sufficient requirement for a non-zero contribution. Being based on observations, the patterns described in eqs.  and 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 \to 1/z$, at symbol level, i.e.  as far as leading-weight terms are concerned. Our basis of respects neither this nor the invariance under complex conjugation at function level: in general ${\mathcal{L}}_\sigma(z,{\bar{z}}) 
\neq {\mathcal{L}}_\sigma(1/z,1/{\bar{z}})$ and ${\mathcal{L}}_\sigma(z,{\bar{z}}) 
\neq {\mathcal{L}}_\sigma({\bar{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 under a standard set of variable transformations. We summarise the most important aspects of these relations in appendix \[app:svhplvariables\]. Quintessentially, these relations determine the coefficients $c_w$ in ${\mathcal{L}}_\sigma(g(z),g({\bar{z}})) = \sum_w c_w {\mathcal{L}}_w(z,{\bar{z}})$ where the sum runs over all words up to weight $|\sigma|$ and, in the present case, $g(x) = 1/x$ or $g(x) = \bar x$.

Let us define \[eq:Fdef\] \_(z,[|[z]{}]{}) 14 ( \_(z,[|[z]{}]{}) + \_([|[z]{}]{},z) + \_(1/z,1/[|[z]{}]{}) + \_(1/[|[z]{}]{},1/z) ) with $\sigma$ 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 ${\mathcal{F}}$s does not span the space of but it does cover the entire space of wavefunctions.

Due to the symmetries of the wavefunction \[eq:wfsyms1\] [\_]{}(z,[|[z]{}]{}) = [\_]{}([|[z]{}]{},z) = [\_]{}(1/z,1/[|[z]{}]{}) = [\_]{}(1/[|[z]{}]{},1/z) and thus \[eq:wfsyms2\] [\_]{}(z,[|[z]{}]{}) = 14 ( [\_]{}(z,[|[z]{}]{}) + [\_]{}([|[z]{}]{},z) + [\_]{}(1/z,1/[|[z]{}]{}) + [\_]{}(1/[|[z]{}]{},1/z) ) one can simply replace ${\mathcal{L}}_\sigma(z,{\bar{z}}) 
\to {\mathcal{F}}_\sigma(z,{\bar{z}})$ to go from the ${\mathcal{L}}$ to the ${\mathcal{F}}$ basis. It may therefore not be immediately obvious how eq.  simplifies the results. Indeed, it requires a few more steps to showcase the advantages of a symmetrised basis.

Firstly, the wavefunction in the ${\mathcal{L}}$ basis contains functions whose indices feature an odd number of the letter $a$. Their leading-weight components are antisymmetric under $z \to 1/z$ because z = - 1/z Converted to ${\mathcal{F}}$ functions they are hence zero at symbol level or, in other words, equal to products of lower-weight 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 ${\mathcal{L}}_\sigma(z,{\bar{z}}) \to {\mathcal{F}}_\sigma(z,{\bar{z}})$

2.  Choose an ${\mathcal{F}}_\sigma(z,{\bar{z}})$ where $\sigma$ contains an odd number of $a$ letters. Plug in definition and rewrite as functions of $z,{\bar{z}}$ using the rules in appendix \[app:svhplvariables\]. The resulting will be of lower weight than the original ${\mathcal{F}}_\sigma$, multiplied by zeta numbers.

3.  Replace again ${\mathcal{L}}_\sigma(z,{\bar{z}}) 
        \to {\mathcal{F}}_\sigma(z,{\bar{z}})$

4.  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 ${\mathcal{L}}_\sigma(z,{\bar{z}}) 
\to {\mathcal{F}}_\sigma(z,{\bar{z}})$ in the wavefunction, *cf.* eqs.  and . To give a few examples for odd-$a$ functions,

$$\begin{aligned}
        {\mathcal{F}}_{a}(z,{\bar{z}}) &= 0 \\
        {\mathcal{F}}_{a,s}(z,{\bar{z}}) = {\mathcal{F}}_{s,a}(z,{\bar{z}}) &= 0 \\
        {\mathcal{F}}_{a,s,s}(z,{\bar{z}}) = {\mathcal{F}}_{s,s,a}(z,{\bar{z}}) &= 4 \zeta_3 \\
        {\mathcal{F}}_{s,a,s}(z,{\bar{z}}) &= -8 \zeta_3 \\
        {\mathcal{F}}_{s,s,s,a}(z,{\bar{z}}) = {\mathcal{F}}_{s,a,s,s}(z,{\bar{z}}) &= 4 \zeta_3 {\mathcal{F}}_s(z,{\bar{z}}).
    \end{aligned}$$

Secondly, we may combine ${\mathcal{F}}_\sigma(z,{\bar{z}})$ and ${\mathcal{F}}_{\tilde \sigma}(z,{\bar{z}})$ with $\tilde \sigma$ the word $\sigma$ reversed, at the cost of generating subleading terms. This is due to the following identity of : \_(z,[|[z]{}]{}) = \_([|[z]{}]{},z) + \_ For a function ${\mathcal{F}}_\sigma$ this entails \_(z,[|[z]{}]{}) = \_(z,[|[z]{}]{}) + \_ \[eq:Fwordrev\] due to the invariance under complex conjugation. Besides removing nearly half of the ${\mathcal{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. – we find

$$\begin{aligned}
    {{{\Omega}^{({1})}}_{\mathrm{2d}}} &= \frac12 C_2 {\mathcal{F}}_s \label{eq:wtd1f} \\
        {{{\Omega}^{({2})}}_{\mathrm{2d}}} &= \frac{1}{8} C_1 C_2 \left({\mathcal{F}}_{a,a}
    -{\mathcal{F}}_{s,s}\right)+\frac{1}{2} C_2^2 {\mathcal{F}}_{s,s} \\
        {{{\Omega}^{({3})}}_{\mathrm{2d}}} &= \frac{1}{16} C_1 C_2^2 \left({\mathcal{F}}_{a,s,a}
    +6 {\mathcal{F}}_{s,a,a}-7 {\mathcal{F}}_{s,s,s}+8 \zeta_3\right) 
    +\frac{1}{16} C_1^2 C_2 
    \left({\mathcal{F}}_{s,s,s}-{\mathcal{F}}_{s,a,a}\right)
    +\frac{3}{4} C_2^3 {\mathcal{F}}_{s,s,s} \\
    \begin{split}
    {{{\Omega}^{({4})}}_{\mathrm{2d}}} &= \frac{1}{16} C_1 C_2^3 \left({\mathcal{F}}_{a,s,s,a}
    +6 {\mathcal{F}}_{s,a,a,s}+4 {\mathcal{F}}_{s,a,s,a}+12 {\mathcal{F}}_{s,s,a,a}
    -23 {\mathcal{F}}_{s,s,s,s}\right.  \\
    &\hspace{4mm} \left.+20 \zeta_3 {\mathcal{F}}_s\right)
    +\frac{1}{64} C_1^2 C_2^2 \left(-{\mathcal{F}}_{a,s,s,a}
    -9 {\mathcal{F}}_{s,a,a,s}-2 {\mathcal{F}}_{s,a,s,a}-24 {\mathcal{F}}_{s,s,a,a}\right.  \\
    &\hspace{4mm} \left.+7 {\mathcal{F}}_{a,a,a,a}+29 {\mathcal{F}}_{s,s,s,s}
    -4 \zeta_3 {\mathcal{F}}_s\right)+\frac{1}{64} C_1^3 C_2 \left({\mathcal{F}}_{s,a,a,s}
    +3 {\mathcal{F}}_{s,s,a,a}\right.  \\
    &\hspace{4mm} \left.-{\mathcal{F}}_{a,a,a,a}-3 {\mathcal{F}}_{s,s,s,s}\right)
    +\frac{3}{2} C_2^4 {\mathcal{F}}_{s,s,s,s} \label{eq:wtd4f}
\end{split}\end{aligned}$$

where we used eq.  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. – to the wavefunction in terms of standard (and the standard $\{0,1\}$ alphabet) in eqs. – shows the benefits of the new basis. In terms of ${\mathcal{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 $-\frac{3}{16} {\mathcal{L}}_1 \zeta_3$ in the coefficient of $C_1^3 C_2$ at four loops .

Finite corrections to the amplitude from two-dimensional evolution {#amplitude}
==================================================================

We now have an algorithm for the calculation of the wavefunction $\Omega_{\rm 2d}$ to any loop order, and we shall use it for the computation of the amplitude. Let us recall from section \[soft\] 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. ]{}. This, in turn, requires the hard part of the two-dimensional wavefunction, which according to [eq. ]{} is obtained by subtracting the $d=2$ limit of the soft wavefunction from the full (two-dimensional) wavefunction $\Omega_{\rm 2d}$ of the previous section. To this end we first define $$\label{Wsoft2d}
    {{\Omega}_{\mathrm{s}}}^{({\rm 2d})}(z,{\bar{z}})\,\equiv \, \left.\lim_{{\epsilon}\to 0} 
    {{\Omega}_{\mathrm{s}}}(p,k)\right|_{\log {\left(\frac{k^2(p-k)^2}{(p^2)^2}\right)}\,\to\, {\mathcal{L}}_s(z,{\bar{z}})}\,,$$ where taking the limit simply corresponds to selecting the leading ${\mathcal{O}}({\epsilon}^0)$ terms in ${{\Omega}_{\mathrm{s}}}(p,k)$. Notice that within the $d=2$ limit we switch to the two-dimensional variables $z$ and ${\bar{z}}$ of eq. , and the single-valued logarithm ${\mathcal{L}}_s(z,{\bar{z}})=\log\frac{z{\bar{z}}}{(1-z)^2(1-{\bar{z}})^2}$ defined in eq. (\[ElDef\]). Having used the symmetrised soft wavefunction we land directly in the class of SVHPLs used in the two-dimensional computation of section \[sec:diffeq\], and in this way the computations of $\Omega_{\rm 2d}$ and ${{\Omega}_{\mathrm{s}}}^{({\rm 2d})}$ are entirely compatible. We note that with the replacement in (\[Wsoft2d\]) the two-dimensional soft wavefunction in  becomes a polynomial in ${\mathcal{L}}_s(z,{\bar{z}})$ at any given order. According to eq. (\[eq:wffullsoftresummed\]), these terms exponentiate and can be resummed to all-orders. Upon applying the changed of variable of (\[Wsoft2d\]) this resummed expression is \[eq:wffullsoftresummed\_SVHPL\] [\_]{}\^[([2d]{})]{}(z,[|[z]{}]{})= \^ e\^[ [C\_2]{}\_s(z,[|[z]{}]{})]{} , with $x = L \, {\alpha_s}/\pi$.

Using the results in section \[2d-bfkl\] for $\Omega_{\rm 2d}$ and the expansion of [eq. ]{} for ${{\Omega}_{\mathrm{s}}}^{({\rm 2d})}$ we determine ${{\Omega}_{\mathrm{h}}}^{({\rm 2d})} = \Omega_{\rm 2d} 
-{{\Omega}_{\mathrm{s}}}^{({\rm 2d})}$, and we can proceed to determine the hard part of the amplitude order by order, according to [eq. ]{}. To this end, recall that the hard wavefunction ${{\Omega}_{\mathrm{h}}}^{({\rm 2d})}$ is guaranteed to integrate to finite terms only, hence it can be integrated in strictly two dimensions. Applying the limit ${\epsilon}\to 0$ in [eq. ]{} to the integrand and the integration measure using the variables $z$ and $\bar{z}$ (*cf.* eqs. (\[k2zzb\]) and (\[eq:dw2d\])) we obtain: \^[(+)]{}\_[NLL,h]{}() = -i [\_[s-u]{}\^2]{}[\^]{}\_[ijij]{}, \[eq:hardAmp2d\] where, in practice, we loop-expand both the wavefunction and amplitude, as was done in eq. . The next two subsections are dedicated to the computation of the integral in eq. (\[eq:hardAmp2d\]), thus determining the hard component of the reduced amplitude order by order. In section \[sec:finiteamp\] we combine the soft and hard components of the reduced amplitude according to eq. (\[eq:redampSplit\]), and finally in section \[sec:hardAmpl\] we similarly combine the soft and hard components of the infrared-renormalized amplitude using eqs. (\[Hsplit\]) and (\[getH2both\]).

To set up the computation of eq. (\[eq:hardAmp2d\]) let us define I [\_]{}\^[([2d]{})]{}(z,[|[z]{}]{}) \[eq:Idef\] 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 $({\rm 2d})$ superscript, and refer to the integrand in (\[eq:Idef\]) as ${{\Omega}_{\mathrm{h}}}(z,{\bar{z}})$. Similarly, while (\[eq:Idef\]) is applied order-by-order, in describing the methods we refrain from using an index for the loop order on either side of (\[eq:Idef\]).

The first method, described in section \[sec:method1\] 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. [@Schwartz:2014wha]. The second method, presented in section \[sec:method2\] below, relies on the symmetry of the wavefunction under inversion, $z \rightarrow 1/z$, ${\bar{z}}\rightarrow 1/{\bar{z}}$, and the action of ${{\hat{H}}_{\mathrm{2d,i}}}$ at fixed external points.

Method I: final integration using the discontinuity of the wavefunction {#sec:method1}
-----------------------------------------------------------------------

Let us define a regularised version of the integral $I$ in eq. : $$I_{\mathrm{reg}}= \frac{1}{4\pi} \int_{\delta^2 < z {\bar{z}}< 1/\delta^2} 
    \frac{{\mathrm{d}}^2 z}{z {\bar{z}}} {{\Omega}_{\mathrm{h}}}(z,{\bar{z}}) \,,
    \label{eq:Iregdef}$$ where the cutoff $\delta$ is assumed to be small, $\delta \ll 1$. The introduction of $\delta$ may seem superfluous at this point as $\lim_{z,{\bar{z}}\to 0} {{\Omega}_{\mathrm{2d}}}= 
\lim_{z,{\bar{z}}\to \infty} {{\Omega}_{\mathrm{2d}}}= \lim_{{\epsilon}\to 0} {{\Omega}_{\mathrm{s}}}$ and thus, using eq. (\[eq:wffullsoftresummed\_SVHPL\]), $\lim_{z,{\bar{z}}\to 0} {{\Omega}_{\mathrm{h}}}= \lim_{z,{\bar{z}}\to \infty} {{\Omega}_{\mathrm{h}}}= 0$; more precisely, ${{\Omega}_{\mathrm{h}}}$ vanishes linearly in $z\bar{z}$ in the soft limit, up to logarithms, rendering the integral in (\[eq:Idef\]) convergent, and the difference $I - I_{\mathrm{reg}}= {\mathcal{O}}(\delta^2)$ (up to logarithms). The necessity of this cutoff despite this good convergence will become clear shortly.

The exclusion of the points $\{0,\infty\}$ in (\[eq:Iregdef\]) enables us to introduce polar coordinates such that $z{\bar{z}}=r^2$ and $\frac{z}{{\bar{z}}}=e^{2i\theta}$, as now all points in the integration region have a non-vanishing Jacobian: $$\label{eq:lastintpol0} I_{\mathrm{reg}}= \frac{1}{4\pi} 
    \int_\delta^{1/\delta} \frac{{\mathrm{d}}r}{r} \int_0^{2\pi} d\theta\,  
    {{\Omega}_{\mathrm{h}}}\left( re^{i\theta},re^{-i\theta} \right)\,.$$ To proceed we express the angular integral in the latter as an integration in the complex $y$ plane where $y\equiv e^{i\theta}$, getting $$\label{eq:lastintpol} I_{\mathrm{reg}}= \frac{1}{4\pi i} \int_\delta^{1/\delta} \frac{{\mathrm{d}}r}{r} 
    \oint_{|y|=1} \frac{{\mathrm{d}}y}{y} \, {{\Omega}_{\mathrm{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 , 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.  runs along the unit circle. Outside of this region, i.e. upon deformation the contour, the in ${{\Omega}_{\mathrm{h}}}(z,{\bar{z}})$ exhibit branch cuts where $z \in [1,\infty]$ and ${\bar{z}}\in [1,\infty]$. In the $r,y$ coordinates of eq.  they correspond to cuts along the real axis in the complex $y$ plane where $y \in [1/r,\infty]$ and $y \in [0,r]$, respectively.

![Position of the branch cuts in $z$ and ${\bar{z}}$ in the complex $y$-plane for $r<1$ (left). The contour along the unit circle in eq.  can be deformed and, consequently, identified with the integral of the ${\bar{z}}$-discontinuity (right), as written in eq. .[]{data-label="fig:ycontour"}](./img/ycontour-crop.pdf)

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 \[fig:ycontour\]. The contour can consequently be shrunk until it corresponds to integrating the ${\bar{z}}$-discontinuity of the wavefunction over $y$ from $0$ to $r$, *cf.* the r.h.s. of figure \[fig:ycontour\]. We can now understand why it is necessary to work with the regularised integral $I_{\mathrm{reg}}$ of eq.  instead of the original $I$ of eq. : while the hard wavefunction ${{\Omega}_{\mathrm{h}}}(z,{\bar{z}})$ vanishes at 0 and $\infty$, its discontinuity, in general, does not. In other words, the contour deformation introduces spurious divergent terms and the cutoff introduced in eq.  regularises them. .

For $r>1$ the branch cuts of $z$ and ${\bar{z}}$ overlap. However, the discontinuity cancels identically in the interval $(1/r,r)$. Repeating the procedure, we again identify the contour integration with integrating the ${\bar{z}}$-discontinuity of ${{\Omega}_{\mathrm{h}}}(z,{\bar{z}})$ over $y$, this time, from $0$ to $1/r$. In total, having modified the contour in (\[eq:lastintpol\]) we find $$\begin{aligned}
    I_{\mathrm{reg}}&= \frac{1}{4\pi i} \left( \int_\delta^1 \frac{{\mathrm{d}}r}{r} \int_0^r \frac{{\mathrm{d}}y}{y} {\mathrm{disc_{{\bar{z}}}}} [{{\Omega}_{\mathrm{h}}}(ry,r/y)] 
    + \int_1^{1/\delta} \frac{{\mathrm{d}}r}{r} \int_0^{1/r} \frac{{\mathrm{d}}y}{y} {\mathrm{disc_{{\bar{z}}}}} [{{\Omega}_{\mathrm{h}}}(ry,r/y)] \right) \nonumber \\
    &= \frac{1}{4\pi i} \left( \int_\delta^1 \frac{{\mathrm{d}}r}{r} \int_1^\infty \frac{{\mathrm{d}}{\bar{z}}}{{\bar{z}}} {\mathrm{disc_{{\bar{z}}}}} [{{\Omega}_{\mathrm{h}}}(r^2/{\bar{z}},{\bar{z}})] 
    + \int_1^{1/\delta} \frac{{\mathrm{d}}r}{r} \int_0^1 \frac{{\mathrm{d}}z}{z} {\mathrm{disc_{{\bar{z}}}}} [{{\Omega}_{\mathrm{h}}}(z,r^2/z)] \right) \nonumber \\
    &= \frac{1}{8\pi i} \left( \int_1^\infty \frac{{\mathrm{d}}{\bar{z}}}{{\bar{z}}} \int_{\delta^2/{\bar{z}}}^{1/{\bar{z}}} \frac{{\mathrm{d}}z}{z} {\mathrm{disc_{{\bar{z}}}}} [{{\Omega}_{\mathrm{h}}}(z,{\bar{z}})] 
    + \int_0^1 \frac{{\mathrm{d}}z}{z} \int_{1/z}^{1/(\delta^2 z)} \frac{{\mathrm{d}}{\bar{z}}}{{\bar{z}}} {\mathrm{disc_{{\bar{z}}}}} [{{\Omega}_{\mathrm{h}}}(z,{\bar{z}})] \right) \nonumber \\
    &= \frac{1}{8\pi i} \left( \int_0^1 \frac{{\mathrm{d}}x}{x} \int_{\delta^2 x}^x \frac{{\mathrm{d}}z}{z} {\mathrm{disc_{{\bar{z}}}}} [{{\Omega}_{\mathrm{h}}}(z,1/x)] 
    + \int_0^1 \frac{{\mathrm{d}}z}{z} \int_{\delta^2 z}^z \frac{{\mathrm{d}}x}{x} {\mathrm{disc_{{\bar{z}}}}} [{{\Omega}_{\mathrm{h}}}(z,1/x)] \right) \label{eq:intdisc}\end{aligned}$$ 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 (\[eq:Iregdef\]) at $z\bar{z}=1$, which admits $(z,{\bar{z}})\leftrightarrow (1/z,1/{\bar{z}})$ symmetry. In the second line of (\[eq:intdisc\]) we reverted to the variable ${\bar{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/{\bar{z}}$ in the first integral and ${\bar{z}}=r^2/z$ in the second; finally in the last line we defined $x=1/{\bar{z}}$ in the both integrals.

![Illustration of the integrations in the $r<1$ (I) and $r>1$ (II) contribution to $I_{\mathrm{reg}}$ of eq.  (white triangles, l.h.s.). They can be viewed as the integral over a square (A) plus two wedges (B) minus two small triangles (C) (delimited by dashed lines, r.h.s.).[]{data-label="fig:method1"}](./img/method1-crop.pdf)

Let us now discuss the evaluation of the final expression in eq. (\[eq:intdisc\]), where the integration region of two terms is depicted as the white area in figure \[fig:method1\]. In order to perform the integration it is useful to view the integrals (*cf.* the r.h.s. of figure \[fig:method1\]) as the integral over a square $$\label{eq:square} I_A(\delta) = \frac{1}{8\pi i} \int_{\delta^2}^1 
    \frac{{\mathrm{d}}z}{z} \int_{\delta^2}^1 \frac{{\mathrm{d}}x}{x} {\mathrm{disc_{{\bar{z}}}}} [{{\Omega}_{\mathrm{h}}}(z,1/x)]\,,$$ plus (the integrals over) two wedges \[eq:wedges\] I\_B() = ( \_0\^1 \_[\^2 x]{}\^[\^2]{} \[[\_]{}(z,1/x)\] + \_0\^1 \_[\^2 z]{}\^[\^2]{} \[[\_]{}(z,1/x)\] ), minus two small triangles \[eq:kite\] I\_C() = ( \_0\^[\^2]{} \_[\^2 x]{}\^x \[[\_]{}(z,1/x)\] + \_0\^[\^2]{} \_[\^2 z]{}\^z \[[\_]{}(z,1/x)\] ), where both $z$ and $x$ are small. Next we would like to evaluate each of these contributions, distinguishing between finite, $\delta$-independent terms, and logarithmically divergent cut-off dependent terms.

The discontinuity w.r.t. ${\bar{z}}$ of ${{\Omega}_{\mathrm{h}}}(z,1/x)$ evaluates to of $z$ and $x$. $I_A(\delta)$ of eq.  thus immediately evaluates to at 1, giving rise to , and at $\delta^2$; the latter contain logarithmically divergent terms in $\delta$. The first (second) integral in the expression of $I_B(\delta)$ in  is calculated close to $z=0$ ($x=0$), *cf.*  figure \[fig:method1\]. 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 \delta^2$, $\log \delta^2 x = 
\log x + \log \delta^2$ and $\log \delta^2 z = \log z + \log \delta^2$, respectively. The outer integrals thereupon generate from their upper limits; in addition it contains logarithmically divergent terms in $\delta$. Contributions from the lower integration limits are dropped according to the (standard) regularisation of : \_[z 0]{} z = 0. \[eq:hplreg\] A similar analysis of $I_C(\delta)$ in eq.  reveals that only powers of $\log \delta^2$ are generated by the integrations over the two small triangles in figure \[fig:method1\].

Since the original integral $I$ of eq. (\[eq:Idef\]) is finite and $I_{\mathrm{reg}}\to I$ for $\delta \to 0$ all terms proportional to $\log \delta^2$ have to cancel between the three contributions $I_A(\delta)$, $I_B(\delta)$ and $I_C(\delta)$. This enables us to derive a simplified integral in which the logarithmically divergent terms are absent altogether whilst giving the same finite terms: $$\begin{gathered}
    I_{\mathrm{reg}}= \frac{1}{8\pi i} \int_0^1 \frac{{\mathrm{d}}z}{z} 
    \int_0^1 \frac{{\mathrm{d}}x}{x} {\mathrm{disc_{{\bar{z}}}}} [{{\Omega}_{\mathrm{h}}}(z,1/x)] \\
    + \frac{1}{8\pi i} \left( - \int_0^1 \frac{{\mathrm{d}}x}{x} 
    \int_0^x \frac{{\mathrm{d}}z}{z} {\mathrm{disc_{{\bar{z}}=1}}} [{{\Omega}_{\mathrm{h}}}(z,1/x)] \big|_{z \ll 1} \right. \\
    \left. - \int_0^1 \frac{{\mathrm{d}}z}{z} \int_0^z \frac{{\mathrm{d}}x}{x} 
    {\mathrm{disc_{{\bar{z}}=1}}} [{{\Omega}_{\mathrm{h}}}(z,1/x)] \big|_{x \ll 1} \right)\,, \label{eq:Iregsimplified}\end{gathered}$$ where all integrals are regulated according to eq.  and ${\mathrm{disc_{{\bar{z}}=1}}} [{{\Omega}_{\mathrm{h}}}(z,1/x)] 
\big|_{z \ll 1}$ and ${\mathrm{disc_{{\bar{z}}=1}}} [{{\Omega}_{\mathrm{h}}}(z,1/x)] \big|_{x \ll 1}$ refer to the aforementioned expansion of the integrand around small $z$ and $x$, respectively. The first integral in eq. (\[eq:Iregsimplified\]) reproduces all finite, cut-off independent terms in $I_A(\delta)$ of (\[eq:square\]), while the second and third ones reproduce, respectively, the finite terms in the two integral in $I_B(\delta)$ in (\[eq:wedges\]); finally, given that no cut-off independent terms are produce by $I_C(\delta)$, it has no trace in (\[eq:Iregsimplified\]).

The above calculation is biased towards the discontinuity with respect to ${\bar{z}}$ 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 ) we expand around $z = {\bar{z}}= 0$ and change variables to the polar coordinates introduced above in eq. . 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 \geq 0$ and $b$ are integer powers. Integrating the azimuth over $[0,\, 2\pi]$ then removes all terms that explicitely depend on $y$, i.e.  that have $b \neq 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 {#sec:method2}
------------------------------------------------------------

The previous method, albeit straightforward on paper, is computationally demanding at high loop orders as it requires extensive use of analytic continuations of 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 \[2d-bfkl\].

Consider the action of ${{\hat{H}}_{\mathrm{2d,i}}}$ on the wavefunction ${{\Omega}_{\mathrm{h}}}(1-z,1-{\bar{z}})$ \[eq:method2int1\] [\_]{}[\_]{}(1-z,1-[|[z]{}]{}) = \^2 w K(w,[|[w]{}]{},z,[|[z]{}]{}) and set $z = {\bar{z}}= 1$ under the integral. Using ${{\Omega}_{\mathrm{h}}}(0,0) = 0$ one gets on the right-hand side: \_[z,[|[z]{}]{}1]{} K(w,[|[w]{}]{},z,[|[z]{}]{}) = K(w,[|[w]{}]{},1,1) [\_]{}(1-w,1-[|[w]{}]{}) \[eq:integrandat1\] with the kernel K(w,[|[w]{}]{},1,1) = + - , *cf.* eq. . It thus follows that (\[eq:method2int1\]), taken in the limit $z,{\bar{z}}\to 1$, yields: $$\begin{aligned}
    \label{eq:method2int1.5} 
    \begin{split}
        \hspace*{-10pt}&{{\hat{H}}_{\mathrm{2d,i}}}{{\Omega}_{\mathrm{h}}}(1-z,1-{\bar{z}}) \big|_{z,{\bar{z}}\to 1} =\\ 
        &\hspace*{20pt} =\, \frac{1}{4\pi} \int {\mathrm{d}}^2 w \left[\frac{ {{\Omega}_{\mathrm{h}}}(1-w,1-{\bar{w}}) }{w{\bar{w}}(1-w)(1-{\bar{w}})} 
        + \frac{ {{\Omega}_{\mathrm{h}}}(1-w,1-{\bar{w}}) }{(1-w)(1-{\bar{w}})} - \frac{ {{\Omega}_{\mathrm{h}}}(1-w,1-{\bar{w}})  }{w{\bar{w}}} \right] 
        \\ &\hspace*{20pt} =\, \frac{1}{4\pi} \int \frac{{\mathrm{d}}^2 w}{w {\bar{w}}} 
        \left[ {{\Omega}_{\mathrm{h}}}\left( \frac{1}{1-w},\frac{1}{1-w} \right) + {{\Omega}_{\mathrm{h}}}(w,{\bar{w}}) 
        - {{\Omega}_{\mathrm{h}}}(1-w,1-{\bar{w}}) \right]\,,
    \end{split}\end{aligned}$$ where in the second line we changed the integration variables in the first two terms – in the first using $w \rightarrow w/(w-1)$, and in the second using $w \rightarrow 1-w$, and then factored out a common denominator. Given that the wavefunction is symmetric under inversion, ${{\Omega}_{\mathrm{h}}}(1/w,1/{\bar{w}}) = {{\Omega}_{\mathrm{h}}}(w,{\bar{w}})$, the first and third terms in the last equation cancel and we find \[eq:method2int3\] [\_]{}[\_]{}(1-z,1-[|[z]{}]{}) |\_[z,[|[z]{}]{}1]{} = [\_]{}(w,[|[w]{}]{}) = I, which can be readily identified with the integral in eq.  which we are interested to compute.

We thus conclude that the integral in eq. , 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 \[sec:wf2d\]. In practice one rewrites the wavefunction ${{\Omega}_{\mathrm{h}}}(1-z,1-{\bar{z}})$ in terms of of $z$ and ${\bar{z}}$, then applies the Hamiltonian by solving the corresponding differential equations, and finally evaluates the resulting expression at $z,{\bar{z}}= 1$. The last step produces the anticipated .

Method I, described in section \[sec:method1\], and method II outlined here show perfect agreement when applied to the wavefunction. However, we emphasise that while former may be applied on individual , the latter can only be applied to expressions which are symmetric under inversion of their arguments, *cf.* eqs. (\[eq:method2int1.5\]) and (\[eq:method2int3\]).

Results for the reduced amplitude {#sec:finiteamp}
---------------------------------

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 \[chap:bfkl\]. There, we claimed that the $\ell$-loop amplitude (divided by $B_0^\ell$, ) has two important number-theoretic properties: all of its terms have weight $\ell$ and there are no terms proportional to $\zeta_2$. We proved this statement for contributions from the soft limit in section \[soft\], see below eq. (\[eq:MsoftExpanded\_8\]). We now show that it holds also for the hard contributions.

We begin by noting that the integrand in (\[eq:Idef\]) 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 \[sec:method1\] and \[sec:method2\], increase the weight of the functions they act on by one, before evaluating the result at $z = {\bar{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 ${{\hat{H}}_{\mathrm{2d,i}}}$ on the wavefunction after a variable transformation $z \to 1-z$. Changing the variables of an 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.

at $z = {\bar{z}}= 1$ evaluate to multiple zeta values () of the same weight, as discussed following eqs. –. We remind the reader that the $(\ell-1)$-loop wavefunction consists of weight-$(\ell-1)$ and weight-$(\ell-1)$ products of and zeta numbers and conclude that the hard contributions to the $\ell$-loop amplitude therefore have uniform weight $\ell$. The absence of $\zeta_2$ in the hard component ${\hat{\mathcal{M}}}_{{\rm NLL,h}}^{(+,\ell)}$ is readily explained by the fact that 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

\[reduced\_hard\_results\] $$\begin{aligned}
    {\hat{\mathcal{M}}}_{{\rm NLL,h}}^{(+,1)} &= 0,  \label{eq:m1hardfinite} \\
    {\hat{\mathcal{M}}}_{{\rm NLL,h}}^{(+,2)} &= 0, \\
    {\hat{\mathcal{M}}}_{{\rm NLL,h}}^{(+,3)}  &= \frac{i\pi}{3!} 
    \bigg\{ \frac{3\zeta_3}{4} {C_1}{C_2}\bigg\} {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}},  \\
    {\hat{\mathcal{M}}}_{{\rm NLL,h}}^{(+,4)} &= 0, \\
    {\hat{\mathcal{M}}}_{{\rm NLL,h}}^{(+,5)} &= \frac{i \pi}{5!} 
    \bigg\{-\frac{5\zeta_5}{2} {C_1}^2 {C_2}^2   
    +\, \frac{45\zeta_5}{2} {C_1}{C_2}^3  \bigg\} {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}},  \\
    {\hat{\mathcal{M}}}_{{\rm NLL,h}}^{(+,6)} &= \frac{i\pi}{6!} 
    \bigg\{\frac{39\zeta_3^2}{16} {C_1}^3 {C_2}^2 
    - \frac{45\zeta_3^2}{2} {C_1}^2 {C_2}^3 
    +\,  \frac{225\zeta_3^2}{2} {C_1}{C_2}^4 \bigg\} {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}},  \\
    {\hat{\mathcal{M}}}_{{\rm NLL,h}}^{(+,7)} &= \frac{i\pi}{7!} 
    \bigg\{ -\frac{2135\zeta_7}{256} {C_1}^4 {C_2}^2 
    + \frac{30135\zeta_7}{256} {C_1}^3 {C_2}^3  
    -\,  \frac{20111\zeta_7}{32} {C_1}^2 {C_2}^4 \nn  \\
    &\hspace{4.0cm}  
    + \frac{6111\zeta_7}{4} {C_1}{C_2}^5\bigg\} {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}, \\
    {\hat{\mathcal{M}}}_{{\rm NLL,h}}^{(+,8)} &= \frac{i\pi}{8!} 
    \bigg\{ \frac{611\zeta_3 \zeta_5}{32} {C_1}^5 {C_2}^2 
    - \frac{643\zeta_3 \zeta_5}{2} {C_1}^4 {C_2}^3 
    + \frac{8597\zeta_3 \zeta_5}{4} {C_1}^3 {C_2}^4 \nn \\ 
    &\hspace{2.0cm}  
    - 7086\zeta_3 \zeta_5 \, {C_1}^2 {C_2}^5    
    +\, 13230 \zeta_3 \zeta_5 \, {C_1}{C_2}^6\bigg\} {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}},
    \label{eq:m8hardfinite} 
    \end{aligned}$$

where we again used the shorthand notation for the colour factors, $C_1=2{C_A}-{{\mathbf{T}}_t^2}$ and $C_2={C_A}-{{\mathbf{T}}_t^2}$. 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. ]{}, i.e. ${\mathcal{H}}^{(+)}_{\rm NLL,h} = {\cal \hat M}^{(+)}_{\rm NLL,h}$, the result of (\[reduced\_hard\_results\]) is sufficient to compute the full infrared-renormalized amplitude ${\mathcal{H}}^{(+)}_{\rm NLL}$ by combining it with the soft contribution ${\mathcal{H}}^{(+)}_{\rm NLL,s}$ of eqs. (\[getH4\]) and (\[eq:HsoftExpanded\_1\]). This will be done in the section \[sec:hardAmpl\] 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. (\[eq:redampSplit\]) and (\[Msh\]), the expressions for the full reduced amplitude through ${\cal O}(\epsilon^0)$ can be easily obtained order-by-order summing the results for the soft amplitude provided in eqs. (\[eq:MsoftExpanded\_1\])-(\[eq:MsoftExpanded\_8\]) and those for the hard amplitude in eqs. (\[eq:m1hardfinite\])–(\[eq:m8hardfinite\]) above, where the former accounts for all infrared singularities plus some finite terms, and the latter for the remaining finite contributions. We obtain

$$\begin{aligned}
    \label{eq:m1finite} 
    {\hat{\mathcal{M}}}_{\rm NLL}^{(1)} &= 
    i\pi {B_{0}} \bigg\{ \frac{1}{2{\epsilon}} \bigg\} {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}, \\ 
    {\hat{\mathcal{M}}}_{\rm NLL}^{(2)} &= 
    i\pi \frac{{B_{0}}^2}{2} \bigg\{ \frac{{C_2}}{4 {\epsilon}^2} \bigg\} {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}, \\  
    {\hat{\mathcal{M}}}_{\rm NLL}^{(3)} &= 
    i\pi \frac{{B_{0}}^3}{3!} 
    \bigg\{ {C_2}^2 \left( \frac{1}{8 {\epsilon}^3} 
    - \frac{11\zeta_3}{4}\right) \bigg\} {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}, \\  
    {\hat{\mathcal{M}}}_{\rm NLL}^{(4)} &= 
    i\pi \frac{{B_{0}}^4}{4!}
    \bigg\{ {C_1}{C_2}^2 \left(
    -\frac{\zeta_3}{8 {\epsilon}} 
    - \frac{3\zeta_4}{16}\right) 
    + {C_2}^3 \left( \frac{1}{16 {\epsilon}^4} 
    + \frac{\zeta_3}{8{\epsilon}} 
    + \frac{3\zeta_4}{16}\right) \bigg\} {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}, \\  \nn
    {\hat{\mathcal{M}}}_{\rm NLL}^{(5)} &= 
    i\pi \frac{{B_{0}}^5}{5!}
    \bigg\{ {C_1}^2 {C_2}^2 \left(
    - \frac{5\zeta_5}{2} \right) 
    + {C_1}{C_2}^3 
    \left( -\frac{\zeta_3}{16 {\epsilon}^2} 
    - \frac{3\zeta_4}{32 {\epsilon}} 
    + \frac{333\zeta_5}{16} \right) \\  
    &\hspace{1.0cm} 
    +\, {C_2}^4 \left( \frac{1}{32 {\epsilon}^5} 
    + \frac{\zeta_3}{16{\epsilon}^2} + \frac{3\zeta_4}{32 {\epsilon}} 
    - \frac{717\zeta_5}{16} \right) \bigg\} {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}, \\ \nn  
    {\hat{\mathcal{M}}}_{\rm NLL}^{(6)} &= 
    i\pi \frac{{B_{0}}^6}{6!}
    \bigg\{ {C_1}^3 {C_2}^2 \bigg( \frac{39 \zeta_3^2}{16} \bigg) 
    + {C_1}^2 {C_2}^3  \bigg( -\frac{399 \zeta_3^2}{16} \bigg) \\ \nn 
    &\hspace{1.0cm} 
    +\, {C_1}{C_2}^4 \bigg(-\frac{\zeta_3}{32 {\epsilon}^3} 
    - \frac{3\zeta_4}{64 {\epsilon}^2}  
    - \frac{3\zeta_5}{32 {\epsilon}}
    + \frac{2637\zeta_3^2}{32}  
    + \frac{5 \zeta_6}{32} \bigg) \\ 
    &\hspace{1.0cm}
    +\, {C_2}^5 \bigg(\frac{1}{64 {\epsilon}^6} 
    + \frac{\zeta_3}{32 {\epsilon}^3} 
    + \frac{3 \zeta_4}{64 {\epsilon}^2}  
    + \frac{3 \zeta_5}{32 {\epsilon}} 
    - \frac{2879 \zeta_3^2}{32} 
    + \frac{5 \zeta_6}{32} \bigg)
    \bigg\} {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}, \\  \nn
    {\hat{\mathcal{M}}}_{\rm NLL}^{(7)} &= 
    i\pi \frac{{B_{0}}^6}{6!}
    \bigg\{ {C_1}^4 {C_2}^2 \bigg( 
    \frac{2135\zeta_7}{256} \bigg)
    {C_1}^3 {C_2}^3 \bigg(  
    \frac{30135\zeta_7}{256} \bigg)
    + {C_1}^2 {C_2}^4 \bigg( \frac{\zeta_3^2}{32 {\epsilon}} 
    + \frac{3\zeta_3 \zeta_4}{32} 
    -\, \frac{20111\zeta_7}{32}\bigg)
    \\ \nn
    &\hspace{1.0cm}
    + {C_1}{C_2}^5 \bigg(
    - \frac{\zeta_3}{64 {\epsilon}^4} 
    - \frac{3\zeta_4}{128 {\epsilon}^3} 
    - \frac{3\zeta_5}{64 {\epsilon}^2}
    - \frac{3\zeta_3^2}{64 {\epsilon}} 
    - \frac{5 \zeta_6}{64 {\epsilon}} 
    - \frac{9 \zeta_3 \zeta_4}{64} 
    +\, \frac{97047 \zeta_7}{64} \bigg)
    \\ \nn
    &\hspace{1.0cm}
    + {C_2}^6 \bigg( \frac{1}{128 {\epsilon}^7}
    + \frac{\zeta_3}{64 {\epsilon}^4} 
    + \frac{3 \zeta_4}{128 {\epsilon}^3} 
    + \frac{3 \zeta_5}{64 {\epsilon}^2}
    + \frac{\zeta_3^2}{64 {\epsilon}} 
    + \frac{5 \zeta_6}{64 {\epsilon}}   
    \\  &\hspace{1.0cm}
    +\, \frac{3 \zeta_3 \zeta_4}{64} 
    - \frac{90711 \zeta_7}{64} \bigg)
    \bigg\} {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}, \\   \nn
    \label{eq:m8finite}
    {\hat{\mathcal{M}}}_{\rm NLL}^{(8)} &= 
    i\pi \frac{{B_{0}}^8}{8!}
    \bigg\{ {C_1}^5 {C_2}^2 \bigg(  
    \frac{611\zeta_3 \zeta_5}{32} \bigg)
    +{C_1}^4 {C_2}^3 \bigg(  
    - \frac{643\zeta_3 \zeta_5}{2} \bigg)
    +{C_1}^3 {C_2}^4 \bigg(  
    \frac{8597\zeta_3 \zeta_5}{4} \bigg) \\ \nn
    &\hspace{1.0cm}
    +\, {C_1}^2 {C_2}^5 \bigg( 
    \frac{\zeta_3^2}{64 {\epsilon}^2} 
    + \frac{3\zeta_3 \zeta_4}{64 {\epsilon}}
    - \frac{228093 \zeta_3 \zeta_5}{32} 
    + \frac{21\zeta_8}{512} \bigg)
    + {C_1}{C_2}^6 \bigg(
    -\frac{\zeta_3}{128 {\epsilon}^5} \\ \nn
    &\hspace{1.0cm}
    -\, \frac{3\zeta_4}{256 {\epsilon}^4}    
    - \frac{3\zeta_5}{128 {\epsilon}^3} 
    - \frac{3 \zeta_3^2}{128 {\epsilon}^2} 
    - \frac{5 \zeta_6}{128 {\epsilon}^2} 
    - \frac{9 \zeta_3 \zeta_4}{128 {\epsilon}} 
    - \frac{9\zeta_7}{128 {\epsilon}} 
    + \frac{749943 \zeta_3 \zeta_5}{64} \\ \nn
    &\hspace{1.0cm} 
    -\, \frac{189 \zeta_8}{1024} \bigg)
    + {C_2}^7 \bigg( \frac{1}{256 {\epsilon}^8} 
    + \frac{\zeta_3}{128 {\epsilon}^5} 
    + \frac{3\zeta_4}{256 {\epsilon}^4}
    + \frac{3\zeta_5}{128 {\epsilon}^3} 
    + \frac{\zeta_3^2}{128 {\epsilon}^2}
    + \frac{5 \zeta_6}{128 {\epsilon}^2}  \\ 
    &\hspace{1.0cm} 
    + \frac{3 \zeta_3 \zeta_4}{128 {\epsilon}} 
    + \frac{9 \zeta_7}{128 {\epsilon}} 
    - \frac{483837 \zeta_3 \zeta_5}{64} 
    + \frac{147 \zeta_8}{1024} \bigg)
    \bigg\} {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}\, ,\end{aligned}$$

These result reproduce the one- to four-loop results of ref. [@Caron-Huot:2013fea] and of our numerically-determined five-loop result in eq. . 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 \[soft\], \[2d-bfkl\] and \[amplitude\].

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. – applies to the amplitude as well: Either of the two methods presented in section \[sec:method1\] and \[sec:method2\] requires us to evaluate at $z = {\bar{z}}= 1$ and we hence expect the presence of (single-valued) starting from weight eleven. Indeed, the eleven-loop amplitude features a term proportional to $g_{5,3,3}$, defined in eq. :

$$\begin{gathered}
    {{{\hat{{\mathcal{M}}}}_{\mathrm{NLL}}}^{(+,{11})}} \supset 
    \frac{1}{102400}\bigg(-\frac{149}{6720} {C_1}^8 {C_2}^2 
    + \frac{26209}{60480} {C_1}^7 {C_2}^3 
    - \frac{14813}{4320} {C_1}^6 {C_2}^4 
    + \frac{210383}{15120} {C_1}^5 {C_2}^5  \\
    - \frac{7549}{252} {C_1}^4 {C_2}^6
    + \frac{39257}{1260} {C_1}^3 {C_2}^7     - 11 {C_1}^2 {C_2}^8 \bigg) \times g_{5,3,3}.
    \label{eq:m11g533}\end{gathered}$$

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 \[soft\]. At twelve loops the reduced amplitude is again comprised of ordinary zeta numbers $\zeta_n$, as there are no weight 12 single-valued MZVs. Such numbers appear then again in the thirteen loop amplitude: $$\begin{gathered}
    {{{\hat{{\mathcal{M}}}}_{\mathrm{NLL}}}^{(+,{13})}} \supset 
    \frac{1}{2207744000}\bigg( \frac{5367943}{497664}  {C_1}^{10} {C_2}^2
    - \frac{32668315}{124416} {C_1}^9 {C_2}^3 
    + \frac{6876365071}{2488320 }  {C_1}^8 {C_2}^4 \\
    - \frac{10213439791}{622080}  {C_1}^7 {C_2}^5 
    + \frac{37444840199}{622080} {C_1}^6 {C_2}^6
    - \frac{10827306157}{77760}  {C_1}^5 {C_2}^7 \\
    + \frac{3841520891}{19440}  {C_1}^4 {C_2}^8
    - \frac{503783639}{3240}  {C_1}^3 {C_2}^9
    + 50459 {C_1}^2 {C_2}^{10} \bigg) \times g_{5,5,3}, 
    \label{eq:m13g553}\end{gathered}$$ and $$\begin{gathered}
    {{{\hat{{\mathcal{M}}}}_{\mathrm{NLL}}}^{(+,{13})}} \supset 
    \frac{1}{2649292800}\bigg( -\frac{1819475}{82944}  {C_1}^{10} {C_2}^2
    + \frac{5621717}{10368}  {C_1}^9 {C_2}^3 
    - \frac{961202489}{165888}  {C_1}^8 {C_2}^4 \\
    + \frac{482408111}{13824}  {C_1}^7 {C_2}^5 
    - \frac{5356152533}{41472}  {C_1}^6 {C_2}^6
    + \frac{1551101681}{5184}  {C_1}^5 {C_2}^7 \\
    - \frac{543921901}{1296} {C_1}^4 {C_2}^8
    + \frac{69045265}{216}  {C_1}^3 {C_2}^9
    - 96967 {C_1}^2 {C_2}^{10} \bigg) \times g_{7,3,3},
    \label{eq:m13g733}\end{gathered}$$ where the single-values zeta numbers $g_{5,5,3}$ and $g_{7,3,3}$ have been defined in eqs.  and .

The fact that the MZV terms in eqs. ,  and  appear in the eleven- and 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 \[soft\].

The infrared-renormalized amplitude\[sec:hardAmpl\]
---------------------------------------------------

We conclude this section by discussing the perhaps most physically relevant infrared-renormalized amplitude (or hard function), which according to [eq. ]{}, is obtained by summing the soft component, given to all orders by the closed expression in [eq. ]{}, and the hard component, which according to [eq. ]{}, coincides with the hard component of the reduced amplitude, ${\cal \hat M}^{(+)}_{\rm NLL,h}$. The latter can be determined to any loop order by following the methods discussed in sections \[sec:method1\] and \[sec:method2\], 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 ${\cal H}$ defined in eq. (\[IRfacteq\]) and loop-expanded following eq. (\[MhatEven\])), up to eight loops:

\[eq:HExpanded\]

$$\begin{aligned}
    \label{eq:HExpanded_1} 
    {\cal H}_{\rm NLL}^{(1)} &= 0, \\ 
    {\cal H}_{\rm NLL}^{(2)} &= 0, \\  
    {\cal H}_{\rm NLL}^{(3)} &= \frac{i\pi}{3!}  
    \bigg\{-{C_2}^2 \frac{11\zeta_3}{4} \bigg\} {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}, \\  
    {\cal H}_{\rm NLL}^{(4)} &= \frac{i\pi}{4!}
    \bigg\{ - C_A {C_2}^2 \frac{3\zeta_4}{16} \bigg\} {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}, \\ 
    {\cal H}_{\rm NLL}^{(5)} &= \frac{i\pi}{5!}
    \bigg\{- C_A^2 {C_2}^2 \frac{5\zeta_5}{2}
    + C_A {C_2}^3 \frac{253\zeta_5}{16}
    -{C_2}^4 \frac{53\zeta_5}{2} \bigg\} {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}, \\ \nn 
    {\cal H}_{\rm NLL}^{(6)} &= \frac{i\pi}{6!}
    \bigg\{ C_A^3 {C_2}^2  \frac{39 \zeta_3^2}{16}
    - C_A^2 {C_2}^3  \frac{141 \zeta_3^2}{8} 
    + C_A {C_2}^4 \bigg(\frac{1275\zeta_3^2}{32} 
    -\, \frac{5\zeta_6}{32} \bigg) 
    \\      & \hspace{2.0cm}
    - {C_2}^5 \frac{481 \zeta_3^2}{16} \bigg\} {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}, \\ \nn 
    {\cal H}_{\rm NLL}^{(7)} &= \frac{i\pi}{7!}
    \bigg\{ - C_A^4 {C_2}^2 \frac{2135 \zeta_7}{256} 
    + C_A^3 {C_2}^3 \frac{21595\zeta_7}{256} 
    + C_A^2 {C_2}^4 \bigg( \frac{3\zeta_3 \zeta_4}{32} 
    -\frac{83293\zeta_7}{256} \bigg) 
    \\
    & \hspace{2.0cm}
    +\, C_A {C_2}^5 \bigg( \frac{3\zeta_3 \zeta_4}{64}
    + \frac{148277 \zeta_7}{256} \bigg)  
    - {C_2}^6 \frac{13443 \zeta_7}{32} \bigg\} {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}, \\   \nn
    \label{eq:HExpanded_8}
    {\cal H}_{\rm NLL}^{(8)} &= \frac{i\pi}{8!}
    \bigg\{ C_A^5 {C_2}^2 \frac{611 \zeta_3 \zeta_5}{32} 
    -C_A^4 {C_2}^3 \frac{7233 \zeta_3 \zeta_5}{32} 
    +C_A^3 {C_2}^4 \frac{16867 \zeta_3 \zeta_5}{16} 
    +C_A^2 {C_2}^5 \bigg(
    - \frac{77383 \zeta_3 \zeta_5}{32} \\
    &\hspace{1.0cm} 
    +\, \frac{21\zeta_8}{512} \bigg)  
    + C_A {C_2}^6 \bigg(
    \frac{174033 \zeta_3 \zeta_5}{64}  
    - \frac{105\zeta_8}{1024} \bigg) 
    - {C_2}^7 \frac{35941\zeta_3 \zeta_5}{32}\bigg\} {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}\, .\end{aligned}$$

Given that ${\mathcal{H}}^{(+)}_{\rm NLL,h} = 
{\cal \hat M}^{(+)}_{\rm NLL,h}$, the same number-theory properties discussed at the end of section \[sec:finiteamp\] apply to the infrared-renormalized amplitude as well. In particular, resummation in terms of gamma functions is excluded.

Numerical analysis and convergence properties {#numerics}
=============================================

The calculation developed in sections \[soft\], \[2d-bfkl\] and \[amplitude\] 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 \equiv L \,{\alpha_s}/\pi$.

Wavefunction
------------

Let us begin by analysing the wavefunction. Given its finiteness, we consider here the leading term in the ${\epsilon}$ expansion, i.e. the two-dimensional soft, hard and full wavefunctions, defined respectively in eqs. (\[Wsoft2d\]), (\[WhardTwod\]) and (\[eq:Htdaction\]).

![Soft, hard and full wavefunction in the complex plane ${\rm Re}(z)$, ${\rm Im}(z)$. Here we plot the singlet component of the wavefunction. The soft and the full wavefunction exhibit singularities at $z = 0$ and $z = \infty$, due to the $z \leftrightarrow 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.[]{data-label="Wave1"}](./img/Omega1-1.pdf "fig:"){width="32.00000%"} ![Soft, hard and full wavefunction in the complex plane ${\rm Re}(z)$, ${\rm Im}(z)$. Here we plot the singlet component of the wavefunction. The soft and the full wavefunction exhibit singularities at $z = 0$ and $z = \infty$, due to the $z \leftrightarrow 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.[]{data-label="Wave1"}](./img/Omega1-2.pdf "fig:"){width="32.00000%"} ![Soft, hard and full wavefunction in the complex plane ${\rm Re}(z)$, ${\rm Im}(z)$. Here we plot the singlet component of the wavefunction. The soft and the full wavefunction exhibit singularities at $z = 0$ and $z = \infty$, due to the $z \leftrightarrow 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.[]{data-label="Wave1"}](./img/Omega1-3.pdf "fig:"){width="32.00000%"}\
![Soft, hard and full wavefunction in the complex plane ${\rm Re}(z)$, ${\rm Im}(z)$. Here we plot the singlet component of the wavefunction. The soft and the full wavefunction exhibit singularities at $z = 0$ and $z = \infty$, due to the $z \leftrightarrow 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.[]{data-label="Wave1"}](./img/Omega1-4.pdf "fig:"){width="32.00000%"} ![Soft, hard and full wavefunction in the complex plane ${\rm Re}(z)$, ${\rm Im}(z)$. Here we plot the singlet component of the wavefunction. The soft and the full wavefunction exhibit singularities at $z = 0$ and $z = \infty$, due to the $z \leftrightarrow 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.[]{data-label="Wave1"}](./img/Omega1-5.pdf "fig:"){width="32.00000%"} ![Soft, hard and full wavefunction in the complex plane ${\rm Re}(z)$, ${\rm Im}(z)$. Here we plot the singlet component of the wavefunction. The soft and the full wavefunction exhibit singularities at $z = 0$ and $z = \infty$, due to the $z \leftrightarrow 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.[]{data-label="Wave1"}](./img/Omega1-6.pdf "fig:"){width="32.00000%"}\
![Soft, hard and full wavefunction in the complex plane ${\rm Re}(z)$, ${\rm Im}(z)$. Here we plot the singlet component of the wavefunction. The soft and the full wavefunction exhibit singularities at $z = 0$ and $z = \infty$, due to the $z \leftrightarrow 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.[]{data-label="Wave1"}](./img/Omega1-7.pdf "fig:"){width="32.00000%"} ![Soft, hard and full wavefunction in the complex plane ${\rm Re}(z)$, ${\rm Im}(z)$. Here we plot the singlet component of the wavefunction. The soft and the full wavefunction exhibit singularities at $z = 0$ and $z = \infty$, due to the $z \leftrightarrow 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.[]{data-label="Wave1"}](./img/Omega1-8.pdf "fig:"){width="32.00000%"} ![Soft, hard and full wavefunction in the complex plane ${\rm Re}(z)$, ${\rm Im}(z)$. Here we plot the singlet component of the wavefunction. The soft and the full wavefunction exhibit singularities at $z = 0$ and $z = \infty$, due to the $z \leftrightarrow 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.[]{data-label="Wave1"}](./img/Omega1-9.pdf "fig:"){width="32.00000%"}

As an example, in figures \[Wave1\] and \[Wave27\] we plot the coefficients $\Omega_{\rm s}^{(\rm 2d)(\ell)}$, $\Omega_{\rm h}^{(\rm 2d)(\ell)}$, $\Omega_{\rm 2d}^{(\ell)}$, 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 $$\begin{aligned}
\label{Tt2num}
\begin{array}{ll}
    \text{singlet}:& {\mathbf T}_{t}^2 \, {\cal M}^{[1]} = 0, \\
    \text{27 representation}:\quad\qquad& 
    {\mathbf T}_{t}^2 \, {\cal M}^{[27]}  
    = 2(N_c+1)\, {\cal M}^{[27]}=8\, {\cal M}^{[27]}.
\end{array}\end{aligned}$$ We plot the wavefunction in the complex $z$ plane, for $\bar z = 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 \leftrightarrow 1/z$ symmetry discussed in sections \[sec:wf2d\] and \[sec:asalphabet\], there is an identical singularity at $z = \infty$ (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. (\[eq:Iregdef\])).

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\gg 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 $\ell$ loop order can be described respectively by the leading logarithms in the two limits, $\sim c_\ell \log^\ell (z \bar z)$ and $\sim c'_\ell \log^\ell\big[1/(1-z)^2(1-\bar z)^2\big]$, where both the magnitude and the sign of the coefficients $c_\ell$ and $c_{\ell}'$ depends on the colour representation considered. Concerning the limit $z = 0$, the asymptotic behaviour is entirely determined by the soft wavefunction, given that $\Omega_{\rm h}^{\rm 2d}(0,0) = 0$. We obtain the coefficients $c_{\ell}$ expanding the soft function in [eq. ]{} (compare eqs. (\[eq:wffullsoftresummed\]) and (\[Well-1-ansatz-sym\])). Taking into account [eq. ]{}, we find \[asymptotic0\] [\^[()]{}]{}(z,|z)|\_[z 0]{} c\_ \^(z |z), c\_= . Given that ${({C_A}-{{\mathbf{T}}_t^2})}= C_A = 3$ for the singlet, while ${({C_A}-{{\mathbf{T}}_t^2})}= -C_A -2 = -5$ for the 27 representation, this explains the sign-oscillating behaviour of the wavefunction near $z = \bar z = 0$ for the singlet, and the constant sign of the 27 representation, which can be seen also in figures \[Wave1\] and \[Wave27\].

Determining the coefficients $c_{\ell}'$ is less trivial, given that near $z = 1$ also $\Omega_{\rm h}^{\rm 2d}$ contributes. An analysis of the asymptotic behaviour up to the 14th order allows us to deduce the pattern and extrapolate. We find: \[asymptotic1\] [\^[()]{}]{}(z,|z)|\_[z 1]{} c’\_ \^ , c’\_= . 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. ]{}. We obtain $$\begin{aligned}
 \label{asymptotic-resum0}
{\Omega}(p,k)|_{z \to 0} &= \frac{{\alpha_s}}{\pi}  \, 
(z \bar z)^{\frac{{\alpha_s}}{2 \pi} \, L \, {({C_A}-{{\mathbf{T}}_t^2})}}, \\
\label{asymptotic-resum1}
{\Omega}(p,k)|_{z \to 1} &=  \frac{{\alpha_s}}{\pi}  \,
\left\{
\begin{array}{ll}
    {}_1F_1\left(2 - \frac{2C_A}{{{\mathbf{T}}_t^2}}, 1, -\frac{{\alpha_s}}{4\pi}\, L \, {{\mathbf{T}}_t^2}\log\left( \frac{1}{(1-z)^2(1-\bar z)^2}\right)\right), & \mbox{if } {{\mathbf{T}}_t^2}\neq 0,  \\
    {}_0F_1\left(1,\frac{{\alpha_s}}{2\pi}\, L \, C_A 
    \log\left( \frac{1}{(1-z)^2(1-\bar z)^2}\right)\right), & \mbox{if } {{\mathbf{T}}_t^2}= 0.
\end{array}
\right.\end{aligned}$$ where ${}_0F_1$ and ${}_1F_1$ are the confluent and Kummer’s confluent hypergeometric function. These resummed expressions are valid only in the leading logarithmic approximation in $z{\bar{z}}$ and $(1-z)(1-{\bar{z}})$, respectively. The generalization of (\[asymptotic-resum0\]) to include subleading logarithms of $z{\bar{z}}$ has been given in (\[eq:wffullsoftresummed\]), while a closed form generalization of (\[asymptotic-resum1\]) is yet unknown.

![Soft, hard and full wavefunction in the complex plane ${\rm Re}(z)$, ${\rm Im}(z)$. Here we plot the component corresponding to the $27$ colour representation. []{data-label="Wave27"}](./img/Omega27-1.pdf "fig:"){width="32.00000%"} ![Soft, hard and full wavefunction in the complex plane ${\rm Re}(z)$, ${\rm Im}(z)$. Here we plot the component corresponding to the $27$ colour representation. []{data-label="Wave27"}](./img/Omega27-2.pdf "fig:"){width="32.00000%"} ![Soft, hard and full wavefunction in the complex plane ${\rm Re}(z)$, ${\rm Im}(z)$. Here we plot the component corresponding to the $27$ colour representation. []{data-label="Wave27"}](./img/Omega27-3.pdf "fig:"){width="32.00000%"}\
![Soft, hard and full wavefunction in the complex plane ${\rm Re}(z)$, ${\rm Im}(z)$. Here we plot the component corresponding to the $27$ colour representation. []{data-label="Wave27"}](./img/Omega27-4.pdf "fig:"){width="32.00000%"} ![Soft, hard and full wavefunction in the complex plane ${\rm Re}(z)$, ${\rm Im}(z)$. Here we plot the component corresponding to the $27$ colour representation. []{data-label="Wave27"}](./img/Omega27-5.pdf "fig:"){width="32.00000%"} ![Soft, hard and full wavefunction in the complex plane ${\rm Re}(z)$, ${\rm Im}(z)$. Here we plot the component corresponding to the $27$ colour representation. []{data-label="Wave27"}](./img/Omega27-6.pdf "fig:"){width="32.00000%"}\
![Soft, hard and full wavefunction in the complex plane ${\rm Re}(z)$, ${\rm Im}(z)$. Here we plot the component corresponding to the $27$ colour representation. []{data-label="Wave27"}](./img/Omega27-7.pdf "fig:"){width="32.00000%"} ![Soft, hard and full wavefunction in the complex plane ${\rm Re}(z)$, ${\rm Im}(z)$. Here we plot the component corresponding to the $27$ colour representation. []{data-label="Wave27"}](./img/Omega27-8.pdf "fig:"){width="32.00000%"} ![Soft, hard and full wavefunction in the complex plane ${\rm Re}(z)$, ${\rm Im}(z)$. Here we plot the component corresponding to the $27$ colour representation. []{data-label="Wave27"}](./img/Omega27-9.pdf "fig:"){width="32.00000%"}

Convergence of the loop expansion of the infrared-renormalized amplitude
------------------------------------------------------------------------

![Partial sums of the *soft* component of the infrared-renormalized amplitude coefficients $\Xi_{\rm NLL,s}^{(+,\ell)}$, up to 15th order, for the singlet (upper plot) and 27 colour representation (lower plot). The dashed vertical line represents the radius of convergence, $R$, determined from the resummed expression.[]{data-label="Radius-Soft"}](./img/Radius1s.pdf "fig:"){width="70.00000%"}\
![Partial sums of the *soft* component of the infrared-renormalized amplitude coefficients $\Xi_{\rm NLL,s}^{(+,\ell)}$, up to 15th order, for the singlet (upper plot) and 27 colour representation (lower plot). The dashed vertical line represents the radius of convergence, $R$, determined from the resummed expression.[]{data-label="Radius-Soft"}](./img/Radius27s.pdf "fig:"){width="70.00000%"}

![Partial sums of the *hard* component of the infrared-renormalized amplitude coefficients $\Xi_{\rm NLL,h}^{(+,\ell)}$, up to 13th order, for the singlet (upper plot) and 27 colour representation (lower plot). The dashed vertical line represents the radius of convergence, $R$, determined by the pole closest to $x = 0$, using the method of Padé approximants.[]{data-label="Radius-Hard"}](./img/Radius1h.pdf "fig:"){width="70.00000%"}\
![Partial sums of the *hard* component of the infrared-renormalized amplitude coefficients $\Xi_{\rm NLL,h}^{(+,\ell)}$, up to 13th order, for the singlet (upper plot) and 27 colour representation (lower plot). The dashed vertical line represents the radius of convergence, $R$, determined by the pole closest to $x = 0$, using the method of Padé approximants.[]{data-label="Radius-Hard"}](./img/Radius27h.pdf "fig:"){width="70.00000%"}

![Partial sums of the infrared-renormalized amplitude coefficients $\Xi_{\rm NLL}^{(+,\ell)}$, up to 13th order, for the singlet (upper plot) and the 27 colour representation (lower plot). The dashed vertical line represents the radius of convergence, $R$, determined by the pole closest to $x = 0$, using the method of Padé approximants.[]{data-label="Radius-Full"}](./img/Radius1.pdf "fig:"){width="70.00000%"}\
![Partial sums of the infrared-renormalized amplitude coefficients $\Xi_{\rm NLL}^{(+,\ell)}$, up to 13th order, for the singlet (upper plot) and the 27 colour representation (lower plot). The dashed vertical line represents the radius of convergence, $R$, determined by the pole closest to $x = 0$, using the method of Padé approximants.[]{data-label="Radius-Full"}](./img/Radius27.pdf "fig:"){width="70.00000%"}

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 fixed-logarithmic accuracy, namely considering the amplitude as a function of $x \equiv L\, {\alpha_s}/\pi$. The high-energy limit adds an interesting twist to the question of convergence, since within the a priori “perturbative regime” where $\alpha_s(\mu^2)$ is small (recall that $\mu^2$ is naturally determined by the momentum transfer $-t$, and we assume $s\gg-t\gg\Lambda_{\rm QCD}^2$) high-energy logarithms can be arbitrarily large, but then the effective expansion parameter $x \equiv L\, {\alpha_s}/\pi$ 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 [@Caron-Huot:2017zfo] 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 ${\cal H}_{\rm NLL}^{(+)}$, which we determined analytically to the 13th order in section \[amplitude\]. For convenience, we introduce the amplitude $\Xi$ and its coefficients $\Xi^{(\ell)}$, defined through $$\begin{aligned}
 \label{MhatEvenB}
    \begin{split}
        {\cal H}_{\rm NLL}^{(+)}  & 
        = \frac{i \pi}{L} \, \Xi_{\rm NLL}^{(+)} \, {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}\,  \\
        &  = \frac{1}{L} \sum_{\ell=1}^\infty  x^\ell  {\cal H}_{\rm NLL}^{(+,\ell)}
        = \frac{i \pi}{L} \sum_{\ell=1}^\infty x^\ell \, \Xi_{\rm NLL}^{(+,\ell)} \, {{\mathbf{T}}_{s-u}^2}{{\mathcal{M}}^{\mathrm{(tree)}}}\, , 
    \end{split}\end{aligned}$$ such that \[MhatEvenC\] [H]{}\_[NLL]{}\^[(+,)]{} = i \_[NLL]{}\^[(+,)]{} [\_[s-u]{}\^2]{}[\^]{} , [eqs.  and ]{}. 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 $\Xi_{\rm NLL}^{(+,\ell)}$ up to thirteen loops can be obtained starting from the analytic expressions given in the ancillary files[^7], using the relations in [eq. ]{}, and converting the multiple zeta values there into decimal numbers. We arrive at \[Xi1\] \_[NLL]{}\^[(+)\[1\]]{} &=& -4.959 x\^3 - 0.2283 x\^4 - 9.230 x\^5 - 2.690 x\^6 - 13.13 x\^7 + 1.696 x\^8\
&&- 20.44 x\^9 + 16.54 x\^[10]{} - 35.99 x\^[11]{} + 46.06 x\^[12]{} - 74.05 x\^[13]{} + (x\^[14]{}),\
\[Xi27\] \_[NLL]{}\^[(+)\[27\]]{} &=& -13.77 x\^3 - 0.6342 x\^4 - 199.2 x\^5 + 381.1 x\^6 - 2826 x\^7 + 9380 x\^8\
&& - 46488 x\^9 + 180393 x\^[10]{} - 797524 x\^[11]{} + 3.239 10\^[6]{} x\^[12]{} - 1.374 10\^[7]{} x\^[13]{}\
&&+ (x\^[14]{}). We consider also the soft and hard contribution to the infrared-renormalized amplitude ${\cal H}_{\rm NLL}^{(+)}$, defined by the two terms in [eq. ]{}. Defining the soft $\Xi_{\rm NLL,s}^{(+)}$ and hard $\Xi_{\rm NLL,h}^{(+)}$, in analogy to [eqs.  and ]{}, we can easily obtain a numerical expression for the singlet and 27 colour representation, as in [eqs.  and ]{}: \[Xis1\] \_[NLL,s]{}\^[(+)\[1\]]{} &=& -7.663 x\^3 - 0.2283 x\^4 - 33.73 x\^5 - 78.04 x\^6 - 210.0 x\^7 - 726.9 x\^8\
&&- 2023 x\^9 - 6237 x\^[10]{} - 18605 x\^[11]{} - 55822 x\^[12]{} - 167566 x\^[13]{} + (x\^[14]{}),\
\[Xis27\] \_[NLL,s]{}\^[(+)\[27\]]{} &=& -15.28 x\^3 - 0.6342 x\^4 - 245.7 x\^5 + 641.8 x\^6 - 4445 x\^7 + 19735 x\^8\
&&- 103863 x\^9 + 507855 x\^[10]{} - 2.566 10\^[6]{} x\^[11]{} + 1.277 10\^[7]{} x\^[12]{}\
&& - 6.398 10\^[7]{} x\^[13]{} + (x\^[14]{}), for the soft part of the infrared-renormalized amplitude, and \[Xih1\] \_[NLL,h]{}\^[(+)\[1\]]{} &=& 2.705 x\^3 + 24.50 x\^5 + 75.34 x\^6 + 196.9 x\^7 + 728.6 x\^8 + 2003 x\^9\
&&+ 6254 x\^[10]{} + 18570 x\^[11]{} + 55869 x\^[12]{} + 167492 x\^[13]{} + (x\^[14]{}),\
\[Xih27\] \_[NLL,h]{}\^[(+)\[27\]]{} &=& 1.503 x\^3 + 46.45 x\^5 - 260.6 x\^6 + 1619 x\^7 - 10356 x\^8 + 57375 x\^9\
&&- 327462 x\^[10]{} + 1.768 10\^[6]{} x\^[11]{} - 9.527 10\^[6]{} x\^[12]{}\
&&+ 5.024 10\^[7]{} x\^[13]{} + (x\^[14]{}), 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 \[Radius-Soft\], \[Radius-Hard\], and \[Radius-Full\].

Considering the all-order resummed expression for the soft component of the infrared-renormalized amplitude in [eq. ]{}, 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. ]{}, and in particular the explicit expression for $\hat \Delta_{\rm NLL}^{(+)}$ in [eq. ]{}, 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}-{{\mathbf{T}}_t^2})}\, 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}-{{\mathbf{T}}_t^2})}$ (the subscript “s” refers to the soft part of the infrared-renormalized amplitude). This corresponds to $R_{s} = 1/3 \simeq 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 $\ell$ for any $x<R_s$, and divergence beyond that point, can indeed be confirmed upon inspecting figure \[Radius-Soft\].

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 $\ell$ 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]: \[padeDef\] \_[NLL]{}\^[(+)]{}|\_[,]{} = . With this definition the Padé approximant has two poles at $x_{\pm} = \big(-b_1 \pm \sqrt{b_1^2 - 4 b_2}\big)/(2 b_2)$. The pole closest to the origin provides a prediction for the radius of convergence of the series: $R = {\rm 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 $\ell$ 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 \[tableXminSoft\]. We see that in both cases the pole closest to $x = 0$ ($x^{[1]}_{\rm s,-} 
= 0.333$ for the singlet and $x^{[27]}_{\rm s,+} =  -0.200$ for the 27 representation) approximates very well the exact radius of convergence, $R_s = 1/{({C_A}-{{\mathbf{T}}_t^2})}$.

  $\ell$   singlet                                                          27                                                                    
  -------- ---------------------------------------------------------------- ------------------------------------------------------------------ -- --
  10       $x^{[1]}_{\rm s,-} = 0.335, \quad x^{[1]}_{\rm s,+} =  -0.703$   $x^{[27]}_{\rm s,-} = 0.472, \quad x^{[27]}_{\rm s,+} =  -0.200$      
  11       $x^{[1]}_{\rm s,-} = 0.333, \quad x^{[1]}_{\rm s,+} =  -1.053$   $x^{[27]}_{\rm s,-} = 0.446, \quad x^{[27]}_{\rm s,+} =  -0.200$      
  12       $x^{[1]}_{\rm s,-} = 0.334, \quad x^{[1]}_{\rm s,+} =  -1.866$   $x^{[27]}_{\rm s,-} = 0.428, \quad x^{[27]}_{\rm s,+} =  -0.200$      
  13       $x^{[1]}_{\rm s,-} = 0.333, \quad x^{[1]}_{\rm s,+} =   3.911$   $x^{[27]}_{\rm s,-} = 0.419, \quad x^{[27]}_{\rm s,+} =  -0.200$      

  : Table summarising the values of $x=\frac{\alpha_s}{\pi}L$ at the poles of the Padé Approximants in eq. (\[padeDef\]), considering the *soft* component of the infrared-renormalized amplitude (indicated by the subscript ${\rm s}$) at orders $\ell=10$ through 13 for the singlet and the 27 representation.[]{data-label="tableXminSoft"}

We thus proceed and apply the same method to the hard component of the infrared-renormalized amplitude. The results are summarised in table \[tableXminHard\].

  $\ell$   singlet                                                          27                                                                     
  -------- ---------------------------------------------------------------- ------------------------------------------------------------------- -- --
  10       $x^{[1]}_{\rm h,-} = 0.333, \quad x^{[1]}_{\rm h,+} =  -0.753$   $x^{[27]}_{\rm h,-} = 0.822, \quad x^{[27]}_{\rm h,+} =  -0.176$       
  11       $x^{[1]}_{\rm h,-} = 0.332, \quad x^{[1]}_{\rm h,+} =  -0.856$   $x^{[27]}_{\rm h,-} = 0.096, \quad x^{[27]}_{\rm h,+} =  -0.179$       
  12       $x^{[1]}_{\rm h,-} = 0.333, \quad x^{[1]}_{\rm h,+} =  -1.258$   $x^{[27]}_{\rm h,-} = -4.392, \quad x^{[27]}_{\rm h,+} =  -0.186$      
  13       $x^{[1]}_{\rm h,-} = 0.333, \quad x^{[1]}_{\rm h,+} =  -1.244$   $x^{[27]}_{\rm h,-} =-0.02, \quad x^{[27]}_{\rm h,+} = -0.185$         

  : Table summarising the values of $x=\frac{\alpha_s}{\pi}L$ at the poles of the Padé Approximants in eq. (\[padeDef\]), considering the *hard* component of the infrared-renormalized amplitude (indicated by the subscript ${\rm h}$) at orders $\ell=10$ through 13 for the singlet and the 27 representation.[]{data-label="tableXminHard"}

We observe that for the singlet there is a highly stable nearest pole at $x^{[1]}_{\rm h,-} = 0.333$. For the 27 representation, in turn, the stable pole at $x^{[27]}_{\rm h,+} \simeq -0.19$ is not always the one closest to the origin, due to the wide fluctuations of $x^{[27]}_{\rm h,-}$. Finally, for the complete infrared-renormalized amplitude we summarise the results in table \[tableXminFull\]. Here we find highly stable results: $x^{[1]}_{+} \simeq  -0.66$ and $x^{[27]}_{+} \simeq -0.24$.

We conclude that Padé approximants based on partial sums of order $\ell=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.

  $\ell$   singlet                                              27                                                        
  -------- ---------------------------------------------------- ------------------------------------------------------ -- --
  10       $x^{[1]}_{-} = 1.092, \quad x^{[1]}_{+} =  -0.624$   $x^{[27]}_{-} = 0.393, \quad x^{[27]}_{+} =  -0.236$      
  11       $x^{[1]}_{-} = 1.266, \quad x^{[1]}_{+} =  -0.666$   $x^{[27]}_{-} = 0.437, \quad x^{[27]}_{+} =  -0.237$      
  12       $x^{[1]}_{-} = 1.311, \quad x^{[1]}_{+} =  -0.661$   $x^{[27]}_{-} = 0.367, \quad x^{[27]}_{+} =  -0.238$      
  13       $x^{[1]}_{-} = 1.466, \quad x^{[1]}_{+} = -0.669$    $x^{[27]}_{-} =0.461, \quad x^{[27]}_{+} = -0.239$        

  : Table summarising the values of $x=\frac{\alpha_s}{\pi}L$ at the poles of the Padé Approximants in eq. (\[padeDef\]), considering the *full amplitude* at orders $\ell=10$ through 13 for the singlet and the 27 representation.[]{data-label="tableXminFull"}

The final results of this analysis are summarised in table \[tableXmin\], 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.

  --------------------------------- ------------------------- ----------------------------------------------------- -- --
  Representation                                      singlet                                                    27    
  Colour factors                             $C_1=6$, $C_2=3$                                     $C_1=-2$,$C_2=-5$    
                                               $1/{C_2}=1/3 $                                        $1/{C_2}=-1/5$    
                                                 $\sim 0.333$                                           $\sim-0.19$    
  Full ${\cal H}_{\rm NLL}^{(+)}$     $-2/{C_2}\simeq -0.666$   $1/\big({C_2}- \frac{3}{8}{C_1}\big) \simeq -0.235$    
  --------------------------------- ------------------------- ----------------------------------------------------- -- --

  : Summary table for the radius of convergence $R$ of the expansion of the infrared-renormalized amplitude in powers of $x=\frac{\alpha_s}{\pi}L$, determined by identifying the pole closest to $x= 0$ using Padé approximants. We use the shorthand notation $C_1 = {(2{C_A}-{{\mathbf{T}}_t^2})}$ and $C_2= {({C_A}-{{\mathbf{T}}_t^2})}$.[]{data-label="tableXmin"}

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 \[Radius-Full\] compared to its soft and hard components in figures \[Radius-Soft\] and \[Radius-Hard\], 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. ]{}, exactly cancels against a similar divergence in the hard component, hence the similar values of $R$ for the soft and hard components in table \[tableXmin\]. 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* ${\cal O}(\epsilon^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 \[Wave1\] and \[Wave27\].

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 \[tableXmin\]. 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\gg-t\gg\Lambda_{\rm QCD}^2$. Given that the expansion parameter, $x=\frac{\alpha_s(-t)}{\pi} 
\log\frac{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.

![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={(2{C_A}-{{\mathbf{T}}_t^2})}$ and $C_2={({C_A}-{{\mathbf{T}}_t^2})}$, based on the Padé approximant analysis for $\ell=11$ and $\ell=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. (\[Rfull\_exact\]).[]{data-label="RR"}](./img/HRR1.pdf "fig:"){width="75.00000%"}\
![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={(2{C_A}-{{\mathbf{T}}_t^2})}$ and $C_2={({C_A}-{{\mathbf{T}}_t^2})}$, based on the Padé approximant analysis for $\ell=11$ and $\ell=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. (\[Rfull\_exact\]).[]{data-label="RR"}](./img/HRR27.pdf "fig:"){width="75.00000%"}

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={(2{C_A}-{{\mathbf{T}}_t^2})}$ and $C_2={({C_A}-{{\mathbf{T}}_t^2})}$[^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 \[RR\] 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 \[RR\] 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: $$\begin{aligned}
\begin{split}
    \hspace*{-30pt}R=\min \left\{x_a,x_b \right\}\qquad
    \text{with}&\qquad x_a = \frac{1}{C_2 - \frac{3}{8}C_1}
    = \frac{1}{{({C_A}-{{\mathbf{T}}_t^2})}- \frac{3}{8}{(2{C_A}-{{\mathbf{T}}_t^2})}}, \\
    &\qquad x_b = - \frac{2}{C_2}=- \frac{2}{{({C_A}-{{\mathbf{T}}_t^2})}}.
    \label{Rfull_exact}
\end{split}\end{aligned}$$ 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_1^{[1]} = 6$, $C_2^{[1]} = 3$. One then obtains from eq. (\[Rfull\_exact\]) a radius of convergence of $R=x^{[1]}_b \simeq -0.6667$, in accordance with the result in table \[tableXmin\]. In turn, $x_a$ gives the pole closest to the origin for the 27 representation, namely for $C_1^{[27]} = -2$, $C_2^{[27]} = -5$, where one obtains from eq. (\[Rfull\_exact\]) $x^{[27]}_a \simeq -0.235$, again, in accordance with table \[tableXmin\].

Conclusions {#conclusion}
===========

In this paper we completed the perturbative calculation of $2\to 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. [@Caron-Huot:2017zfo] where we determined the infrared singularities, we now computed the finite corrections to the hard amplitude ${\cal 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 [@Caron-Huot:2017zfo], 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. [@Caron-Huot:2017zfo] 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. (\[getH2\]) and (\[eq:redampSplit\]), each separated part of the BFKL-motivated reduced amplitude needs only be calculated to order ${\cal O}(\epsilon^0)$, and by carefully recombining them we obtained the renormalized amplitude in eq. (\[eq:HExpanded\]). 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 \[2d-bfkl\]. The result lives inside a very rigid space of functions: the $\ell$-loop wavefunction is a linear combination of weight-$\ell$ single-valued harmonic polylogarithms (SVHPLs) of $z$ and ${\bar{z}}$ 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$ and ${\bar{z}}$.

The hard contribution to the infrared-renormalized amplitude ${\cal 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. (\[ReducedAmpNLLresum2B\]), 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. (\[reduced\_hard\_results\]) and (\[eq:m11g533\])). The presence of multi zeta values – which make their first appearance at weight 11 involving a single-valued version [@Brown:2013gia; @Schnetz:2013hqa] of $\zeta_{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 ${\cal H}$ to high loop order, it is straightforward to study the results numerically. In section \[numerics\] 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,{\bar{z}}\to 0$ and $z,{\bar{z}}\to \infty$, while the third $z,{\bar{z}}\to 1$ corresponds to the limit of large internal momentum. The former are described analytically by the soft wavefunction in eq. (\[eq:wffullsoftresummed\]), 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 incomplete – cancellation 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. (\[asymptotic-resum1\])) 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 ${\cal O}(\epsilon^0)$ infrared-renormalized amplitude has a finite radius of convergence in the expansion parameter $x=L\alpha_s/\pi$. For the soft contribution, where we have a resummed analytic expression, eq. (\[ReducedAmpNLLresum2\]), this radius of convergence can readily be identified as the first pole of a gamma function, generating asymptotic behaviour $\sim (x(C_A-{{\mathbf{T}}_t^2}))^{\ell}$ at high orders, $\ell\to \infty$. 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 \[Radius-Full\] (compare with figures \[Radius-Soft\] and \[Radius-Hard\] 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 ${{\mathbf{T}}_t^2}$, eq. (\[Rfull\_exact\]) 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. [@Larkoski:2016zzc]. It remains for future work to understand the true high-energy (large $x=L\alpha_s/\pi$) behaviour.

Let us conclude with a brief summary of the state-of-the-art knowledge of partonic $2\to 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* (\[Mreal\]) with two-loop corrections to the trajectory $\alpha^{(2)}_g(p^2)$, and suitable impact factors (the former, in particular, was calculated in [@Fadin:1995xg; @Fadin:1996tb; @Fadin:1995km; @Blumlein:1998ib]; it can also be extracted from two-loop calculations of $2\to 2$ scattering amplitudes [@DelDuca:2001gu]). 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. [@Caron-Huot:2017fxr], 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 indeed – using the techniques we developed in the present paper – technically feasible, to compute higher-loop corrections in this tower of logarithms. NNLL corrections in the signature-even sector are in turn simpler and can be deduced from linear BFKL evolution with a NLO kernel [@Fadin:1995xg; @Fadin:1996tb; @Fadin:1995km; @Blumlein:1998ib], 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 ${\cal N}=4$ supersymmetric Yang-Mills theory, see e.g. [@Caron-Huot:2016owq; @Dixon:2016nkn], but also, more recently in the context of the singularity structure of gauge theories including QCD [@Almelid:2017qju]. In both cases, the high-energy limit served as crucial input.

We would like to thank Lorenzo Magnea, Jenni Smillie and Claude Duhr for useful discussions. We would also like to thank Gudrun Heinrich and Stephan Jahn for helpful feedback on the numerical evaluations of integrals with pySecDec. EG’s research is supported by the STFC Consolidated Grant ‘Particle Physics at the Higgs Centre’, JR’s research was supported by the Walter Nimmo and Walter Scott PhD studentship, LV research is supported by Fellini - Fellowship for Innovation at INFN, funded by the European Union’s Horizon 2020 research programme under the Marie Skłodowska-Curie Cofund Action, grant agreement no. 754496. EG and JR would like to thank the McGill Physics department for their kind hospitality. EG, SCH and LV would like to thank the GGI, Florence, for support during the programme ‘Amplitudes in the LHC era’ in Autumn 2018. EG thanks the Simons Foundation for support as a Simons Visiting Scientist in GGI in Autumn 2018 and the CERN theory department for hospitality as a Scientific Associate in 2019. This research is also supported by the National Science and Engineering Council of Canada, the Canada Research Chair program, and the Fonds de Recherche du Québec - Nature et Technologies.

Harmonic polylogarithms {#app:hpls}
=======================

Harmonic polylogarithms () [@Remiddi:1999ew] extend the natural logarithm $\log z$ with $z \in \mathbb{C}$ to nested integrals. Similarly to the well-known polylogarithms ${\mathrm{Li}}_n(z)$ they are defined recursively namely \[eq:hpldef\] H\_[0,]{}(z) = \_0\^z t H\_[1,]{}(z) = \_0\^z t where $\sigma$ is a “word” of any length made from the letters[^10] $\{0,1\}$. The number of indices of a $H(z)$ is called the *weight* of the function. By means of eq.  it corresponds to the number of nested integrals. The recursion is closed by the weight-1 identities \[eq:hpl1def\] H\_[0]{}(z) = z H\_[1]{}(z) = -(1-z). form a shuffle algebra and thus obey shuffle product identities \[eq:hplshuffle\] H\_(z) H\_(z) = \_ H\_(z) where $\rho \shuffle \sigma$ denotes the shuffle of the words $\rho$ and $\sigma$.

The indices of a may be shortened by means of a collapsed notation; one replaces strings of zeros followed by a one according to \_[n ]{},1 n+1 \[eq:collapsedn\] for example $H_{0,1,0,0,1,1}(z) \to 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 as nested sums. One commonly used definition is \[eq:hplseriesdef\] H\_(z) = \_[j=1]{}\^z\^j Z\_j() with \[eq:zsumdef\] Z\_j(a,) = \_[i=2]{}\^j Z\_[i-1]{}() Z\_j(1) = 1/j where we assume the collapsed notation. Note that the aforementioned depth is equal to the number of nested sums.

The Taylor series of , defined by eq. , whose rightmost index is non-zero, is given by eq.  with . Trailing zeros in the indices of a 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) and solves for $H_{\sigma,0}(z)$. This procedure can be applied recursivly until all trailing zeros are removed. Hence, can always be written as a series in $z$ and $\log z$.

For arguments between 0 and 1 yield real values. They show branch cuts on the real axis where $z \in [1,\infty)$ and are thus multi-valued functions.

Single-valued harmonic polylogarithms {#app:svhpls}
=====================================

Single-valued harmonic polylogarithms () [@Brown:2004ugm] are the class of all branch cut-free, single-valued, combinations of . Their construction is somewhat involved and we will only provide a short summary here. Further details can be found in e.g. refs. [@Pennington:2012zj; @Dixon:2012yy; @DelDuca:2013lma].

are functions of a complex variable $z$ and its complex conjugate ${\bar{z}}$. They correspond to the linear combinations of $H_\sigma(z) H_{\sigma'}({\bar{z}})$ that solve \[eq:ddzL\] \_[0,]{}(z,[|[z]{}]{}) = \_[1,]{}(z,[|[z]{}]{}) = and obey the boundary conditions [@Pennington:2012zj] \_(z,[|[z]{}]{}) = 1, \_[0\_n]{}(z,[|[z]{}]{}) = \^n (z [|[z]{}]{}) / n! \_[z 0]{} \_[0\_n]{}(z,[|[z]{}]{}) = 0. 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 $$\begin{aligned}
    y_0 &= x_0 \label{eq:y0x0} \\
    \tilde Z(y_0,y_1) y_1 \tilde Z(y_0,y_1)^{-1} 
    &= Z(x_0,x_1)^{-1} x_1 Z(x_0,x_1) \label{eq:y1x1}\end{aligned}$$ where $Z$ is the so-called Drinfeld associator. It is defined as the generating series Z(x\_0,x\_1) = \_[X\^\*]{} H\_(1) Z(y\_0,y\_1) = \_[Y\^\*]{} H\_[()]{}(1) \[eq:drinfeld\] where the “tilde” operation reverses words and $\phi$ maps $y_i \to x_i$. The values of the at $z = 1$ in the definition are regularised by the shuffle algebra. Eq.  can be solved iteratively for $y_1$.

The can then be extracted from the product of another two generating series \_[X\^\*]{} \_(z,[|[z]{}]{}) = L\_X(z) L\_Y([|[z]{}]{}) where L\_X(z) = \_[X\^\*]{} H\_(z) L\_Y([|[z]{}]{}) = \_[Y\^\*]{} H\_[()]{}([|[z]{}]{}) with “tilde” and $\phi$ defined below eq..

obey the same shuffle product as , namely \[svhplshuffle\] \_(z,[|[z]{}]{}) \_(z,[|[z]{}]{}) = \_ \_(z,[|[z]{}]{}).

Holomorphic part and single-value map {#app:holomorphicpart}
-------------------------------------

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 $\psi(z,{\bar{z}})$ as the limit \[eq:holopartdef\] \^[(h)]{}(z) = (z,0) |\_[[|[z]{}]{}0]{}. For a given linear combination of taking this limit simply amounts to replacing ${\mathcal{L}}_\sigma(z,{\bar{z}}) 
\to H_\sigma(z)$.

The dependence on ${\bar{z}}$ is reconstructed by the single-value map \[eq:sdef\] ( \^[(h)]{}(z) ) = (z,[|[z]{}]{}) which is discussed in detail in refs. [@Brown:2013gia; @DelDuca:2016lad]. Again, we restrict ourselves here to stating the (obvious) replacement rule $H_\sigma(z) \to 
{\mathcal{L}}_\sigma(z,{\bar{z}})$ which generates the corresponding single-valued expression from a linear combination of of $z$. As the action of the Hamiltonian ${{\hat{H}}_{\mathrm{2d,i}}}$ removes constant terms from the wavefunction prior to integration we shall not discuss this aspect in the context of eqs.  and here. The interested reader is referred to the above references.

Variable transformations {#app:svhplvariables}
------------------------

obey relations under certain variable transformations. For the most part they are, in some sense, the *same* relations that apply to due to the single-value map discussed above in appendix \[app:holomorphicpart\]. While the latter are much better documented (for an overview we recommend ref. [@Maitre:2005uu]) we struggled to find a comprehensive list for which motivated this appendix.

In section \[sec:asalphabet\] we transform $z \to 1/z$ and $z \leftrightarrow {\bar{z}}$ to account for the symmetries of the two-dimensional wavefunction. In addition, we consider $z \to 1-z$ in section \[sec:method2\] to facilitate the “last integration”. Let us discuss the latter transformation in detail.

At the level of it is straightforward to find relations under $z \to 1-z$. Effectively, the transformation moves the lower limit of the integral definition from zero to one. Consider the weight-$w$ with argument $1-z$ $$\begin{aligned}
H_{0,a_2,\dots,a_w}(1-z) &= \int_0^{1-z} \frac{{\mathrm{d}}t}{t} H_{a_2,\dots,a_w}(t) \nn \\
&= \int_0^1 \frac{{\mathrm{d}}t}{t} H_{a_2,\dots,a_w}(t) - \int_{1-z}^1 \frac{{\mathrm{d}}t}{t} H_{a_2,\dots,a_w}(t) \nn \\
&= H_{0,a_2,\dots,a_w}(1) - \int_0^z \frac{{\mathrm{d}}t}{1-t} H_{a_2,\dots,a_w}(1-t) \label{eq:hpl0zerotoone}\end{aligned}$$ and $$\begin{aligned}
H_{1,a_2,\dots,a_w}(1-z) &= \int_0^{1-z} \frac{{\mathrm{d}}t}{1-t} H_{a_2,\dots,a_w}(t) \nn \\
&= \int_0^1 \frac{{\mathrm{d}}t}{1-t} H_{a_2,\dots,a_w}(t) - \int_{1-z}^1 \frac{{\mathrm{d}}t}{1-t} H_{a_2,\dots,a_w}(t) \nn \\
&= H_{1,a_2,\dots,a_w}(1) - \int_0^z \frac{{\mathrm{d}}t}{t} H_{a_2,\dots,a_w}(1-t) \label{eq:hpl1zerotoone}\end{aligned}$$ with H\_0(1-z) = -H\_1(z) H\_1(1-z) = -H\_0(z). Since the inside the integrals in eqs.  and are of weight $w-1$ this defines a recursive prescription of how to write any of $1-z$ in terms of of $z$.

By means of the holomorphic part of and the single-value map, see appendix \[app:holomorphicpart\], these relations can be applied to . However, it is also possible to solve the recursion and write the answer directly as a sum. We find \_[a\_1,…,a\_w]{}(1-z,1-[|[z]{}]{}) = \_[j=0]{}\^w (-1)\^j \_[a\_1,…,a\_j]{}(z,[|[z]{}]{}) \_[a\_[j+1]{},…,a\_w]{}(1,1) with the “$\sim$” operation swapping the indices $0 \leftrightarrow 1$.

Similarly, on can derive identities for the transformation $z \to 1/z$, ${\bar{z}}\to 1/{\bar{z}}$. Again, the recursion can be solved and the resulting formula is simply yet slightly awkward to write out. To do so we define $n_0(\sigma)$ ($n_1(\sigma)$) to count the number zeros (ones) in the indices $\sigma$ and $\hat s_{1 \to 0+1}$ to split ${\mathcal{L}}_\sigma(z,{\bar{z}})$ into a sum of $2^{n_1(\sigma)}$ according to the index rule $1 \to 0 + 1$. For example, $$\begin{gathered}
\hat s_{1 \to 0+1} \left[ {\mathcal{L}}_{1,0,0,1,0}(z,{\bar{z}}) \right] = \\
{\mathcal{L}}_{0,0,0,0,0}(z,{\bar{z}}) + {\mathcal{L}}_{0,0,0,1,0}(z,{\bar{z}}) 
+ {\mathcal{L}}_{1,0,0,0,0}(z,{\bar{z}}) + {\mathcal{L}}_{1,0,0,1,0}(z,{\bar{z}})\end{gathered}$$ Then \_[a\_1,…,a\_w]{} ( 1z,1[|[z]{}]{}) = \_[j=0]{}\^w (-1)\^[n\_0(a\_1,…,a\_j)]{} s\_[1 0+1]{} \_[a\_[j+1]{},…,a\_w]{}(,). The values of at $z,{\bar{z}}\to \infty$ are related to the values at $z,{\bar{z}}= 1$ by yet another transformation: $z \to z/(z-1)$. \_[a\_1,…,a\_w]{} ( , ) = (-1)\^[n\_1(a\_1,…,a\_w)]{} s\_[0 0+1]{} with $\hat s_{0 \to 0+1}$ defined like $\hat s_{1 \to 0+1}$ but based on the index rule $0 \to 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 at $z,{\bar{z}}= 1$.

Lastly, let us examine the transformation $z \leftrightarrow {\bar{z}}$ and how to related an ${\mathcal{L}}_\sigma({\bar{z}},z)$ to (a sum of) ${\mathcal{L}}_{\sigma_i'}(z,{\bar{z}})$. The easy yet computationally heavy way is to translate ${\mathcal{L}}_\sigma(z,{\bar{z}})$ to , swap $z \leftrightarrow {\bar{z}}$, extract the holomorphic part by means of eq.  and finally apply $\mathbf{s}$ . For 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 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. . Consider the weight-$n$ ${\mathcal{L}}_\sigma(z,{\bar{z}})$ with $\sigma = \sigma_1,\dots,\sigma_n$ and swap $z \leftrightarrow {\bar{z}}$. Then \_([|[z]{}]{},z) = \_(z,[|[z]{}]{}) + \_[i=4]{}\^[||]{} \_[j=0]{}\^[||-i]{} y1(\_j,…,\_[i+j]{}) \_[\_[(A)]{},\_[(B)]{},1]{} \[eq:ztozb\] where the “tilde” map was defined below eq.  and $y1(\sigma)$ is the coefficient of the product of $x_0$ and $x_1$ corresponding to $\sigma$, e.g. if $\sigma 
= 1,1,0,1,0$ then $y1(\sigma)$ is the coefficient of $x_1 x_1 x_0 x_1 x_0$. The indices (A) in eq.  only appear if $j - 1 \geq 1$ and likewise (B) if $i + j + 1 \leq n$.

[^1]: Note that ultraviolet renormalization is irrelevant for the signature-even amplitude at the logarithmic accuracy considered.

[^2]: This is a direct consequence of the fact that we have removed the factor of the gluon Regge trajectory in defining the reduced amplitude in [eq. ]{}.

[^3]: For the most part of this section we will use the standard letters, 0 and 1. Only in section \[sec:asalphabet\] we introduce a new alphabet to simplify the two-dimensional evolution.

[^4]: use the collapsed notation, *cf.* eq.  in appendix \[app:hpls\].

[^5]: Brown [@Brown:2013gia] refers to it as $\zeta_{\rm sv}(3,5,3)$ while Schnetz [@Schnetz:2013hqa] calls it $g_{335}$.

[^6]: This step requires some amount of creativity but is greatly helped by The On-Line Encyclopedia of Integer Sequences (), `https://oeis.org`.

[^7]: The same result is provided explicitly in the main text in [eqs. –]{}, up to eight loops.

[^8]: There is of course some freedom of in choosing the degrees of the polynomials in the numerator and the denominator. After some experimentation we found that Padé approximants with second-order denominators yield stable predictions for the position of the first pole already at relatively low orders, and hence we use this form as the default choice for the analysis presented here. Qualitatively, the results are the same using different Padé approximants.

[^9]: The hard function is provided in [eqs.  and ]{} in terms of $C_A$ and $C_2$, but we find it more convenient for this analysis to express it as a function of the color operators $C_1$ and $C_2$.

[^10]: The full alphabet of includes the letter $-1$. In the present work however we only encounter integrals corresponding to the letters $0$ and $1$.
