Subleading Regge limit from a soft anomalous dimension

Wilson lines capture important features of scattering amplitudes, for example soft effects relevant for infrared divergences, and the Regge limit. Beyond the leading power approximation, corrections to the eikonal picture have to be taken into account. In this paper, we study such corrections in a model of massive scattering amplitudes in N=4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{N}=4 $$\end{document} super Yang-Mills, in the planar limit, where the mass is generated through a Higgs mechanism. Using known three-loop analytic expressions for the scattering amplitude, we find that the first power suppressed term has a very simple form, equal to a single power law. We propose that its exponent is governed by the anomalous dimension of a Wilson loop with a scalar inserted at the cusp, and we provide perturbative evidence for this proposal. We also analyze other limits of the amplitude and conjecture an exact formula for a total cross-section at high energies.

Future physics analyses at the LHC will require conceptional advances in the theoretical understanding of scattering processes. One new frontier will be higher-loop processes depending on many mass and kinematic scales, e.g. when considering mixed QCD and electroweak processes. While in some cases a numerical approach may be feasible and adequate, it seems clear that conceptual breakthroughs will be driven by new analytic ideas. When dealing with processes depending on many scales, an important question is to understand in which situations systematic expansions can be applicable, and how the latter can be obtained systematically. One particularly interesting and important limit is the eikonal limit, which describes emission of soft radiation. At leading power, the physics is described by correlation functions of Wilson lines. Many recent papers are dedicated to studying power corrections to the eikonal limit [2][3][4].
It is often extremely helpful to have a toy model at hand for developing new ideas. Among various Yang-Mills theories, the N = 4 super Yang-Mills theory stands out as a particularly nice model. It is often dubbed the hydrogen atom of quantum field theory, due to a hidden symmetry that is the generalization of the Laplace-Runge-Lenz symmetry of the hydrogen atom. In order to be able to study massive scattering amplitudes, we give a vacuum expectation value to some of the scalar fields in the model. In this way, we can consider fourparticle amplitudes depending on two Mandelstam variables and a mass. The amplitudes are ultraviolet and infrared finite, so that they can be evaluated directly in four dimensions.
One of many examples where this model brought about conceptual advances for generic quantum field theories is in understanding the structure of Feynman integrals [5], which is closely related to their analytic evaluation [6]. The three-loop planar Feynman integrals needed for the amplitude mentioned above were computed in ref. [1], using a version of the differential equations method with improvements for integrals finite in four dimensions. These formulas provide a fully analytic result for a three-loop four-point amplitude depending on three scales.
One exciting feature of the amplitude is that many of its physically interesting limits are either exactly known, or governed by integrability. This is the case for the Regge limit, which at leading power is controlled by the anomalous dimension of a cusped Wilson loop, a problem that is known to be integrable [7,8]. The amplitude is exactly known in the highenergy limit, and equal to the Bern-Dixon-Smirnov (BDS) ansatz [9][10][11]. Moreover, the lowenergy limit is described by an effective action, and in the forward limit the amplitude can be interpreted as a total cross-section of producing massive W bosons and other particles; an exact formula for that cross-section at high energies will be conjectured below. Finally, it is possible to study threshold effects that are related to bound states of a hydrogen-like system [12].
In this paper, we study in detail the various physical limits mentioned above, focusing in particular on the Regge limit. It is well-known that in the planar limit it is described by a single power law, involving the gluon Regge trajectory, which is a nontrivial function of the internal masses and momentum transfer. Surprisingly, we will observe from the explicit three-loop amplitude that the first subleading power is also controlled by a single power JHEP04(2018)047 law. This is rather remarkable in light of the current understanding of subleading powers in the Regge limit, where, for example, quark loops typically produce double logarithmic corrections (see [13,14]), In contrast, in this theory we find only a single power of logarithm per loop order.
In order to better understand this phenomenon and hopefully initiate a systematic expansion in the Regge limit, we will make use of a special property of this model, wherein the Regge limit can be mapped to a limit of a massless internal leg. This was used before to give an alternative definition of the Regge trajectory as the anomalous dimension of a cusped Wilson loop with a finite angle. This kind of soft limit is at the moment under better theoretical control since it is conceptually closer to the limit of soft external momenta studied in [2][3][4]. We will make a proposal for an independent definition of the first subleading exponent, namely as the anomalous dimension of a Wilson loop with a scalar inserted at the cusp. We test this proposal by explicitly computing this anomalous dimension up to two loops.
The paper is organized as follows. In section 2, we define the model and amplitudes under consideration. We discuss in detail the various physical limits and point out the all-loop structure expected in some of them, using the one-loop result as a pedagogical example. Section 3 summarizes our observations about the Regge limit at next-to-leading power, up to three loops. Then, in section 4 we compute the anomalous dimension of a cusped Wilson loop with a scalar insertion, and test our proposal that its anomalous dimension is equal to the second Regge trajectory appearing at subleading power. We present our conclusions in section 5. The paper contains several appendices with technical details. Appendix A collects our results for a total cross-section, and evidence for its conjectured high-energy limit. Appendix B explains the use of dual conformal symmetry to conveniently parametrize the Regge expansion. Appendix C contains a detailed account of the analytic continuation and differential equation technology needed to derive the various expansions of the three-loop scattering amplitude. In appendix D we discuss the calculation of the Feynman integrals for the two-loop Wilson line calculation.
2 Massive scattering amplitudes in N = 4 super Yang-Mills

Setup and four-particle amplitudes
We consider the N = 4 super Yang-Mills theory in the planar limit. We spontaneously break the SU(N c ) gauge group to SU(N c −4)×SU(4)×U (1). In this way, in addition to the "gluons" of the unbroken SU(N c −4) part of the gauge group, we have massless bosons from the unbroken SU(4) subgroup, a U(1) photon, and massive W bosons from the off-diagonal part. In what follows we will take N c large and discuss the leading term of the amplitudes.
As discussed in ref. [11], this allows us to consider color-ordered amplitudes YȲ → YȲ . 1 For further papers discussing various aspects of massive amplitudes on the Coulomb branch of N = 4 super Yang-Mills, see [16][17][18][19][20][21][22][23][24][25][26][27]. Here the particle Y is one of the off-diagonal generators of the unbroken SU(4) subgroup, lying above the diagonal;Ȳ is then the Hermi-JHEP04(2018)047 forward limit t = 0, total cross-section, conjectured exact formula at high energy lim s→∞ σ(s) Regge limit s → ∞, Wilson line description, integrability high energy limit s, t → ∞, exact formula, infrared/collinear divergences regulated by mass threshold expansion s ∼ 4m 2 , relation to hydrogen-like systems soft limit s, t → 0, effective action description Amplitude A(s, t, m 2 ) Figure 1. The scattering amplitude A(s, t, m 2 ) has various physically interesting and overlapping limits. In many of the latter, exact results are known or conjectured (e.g. high-energy limit), while other limits are known to be governed by integrability. tian conjugate. We choose Y to be a complex scalar within the N = 4 supermultiplet, but since the setup preserves supersymmetry any other helicity choice would give an equivalent result. An important motivation for considering such amplitudes is that they are natural from the AdS/CFT viewpoint [28], and that they have a dual conformal symmetry [11].
At tree-level, the result for the scattering amplitude is the same as in the unbroken theory, Amplitudes with other external states, such as gluons, are related to this one by supersymmetry. At loop level, at leading order in N c , the interactions are mediated by massive W bosons, whose mass provides a natural infrared regularization. We define The amplitude M can be expanded perturbatively in the coupling g 2 ≡ g 2 YM N c /(16π 2 ), as

JHEP04(2018)047
The expression for the loop integrand of M up to four loops was derived using unitarity cuts [29]. The loop integrals up to three loops were evaluated analytically in ref. [1]. The main focus of this paper is to investigate the various limits discussed above, and to understand surprising structures appearing in them. We use the technology of ref. [1] to derive the expansions.
In this section, we use the one-loop expressions as a pedagogical example, and point out the all-loop structure, whenever the latter is known. The one-loop term M (1) of eq. (2.3) is given by a massive one-loop box integral, which evaluates to (the form below is due to [30]) Here we introduced dimensionless variables 5) and the following abbreviations, The functions appearing in eq. (2.4) are examples of polylogarithms, with the dilogarithm defined as Li 2 (x) = − x 0 log(1 − y)/y dy. The above formulas are valid in the Euclidean region u, v > 0. In order to continue to other regions, a small imaginary part has to be added to s and t, according to the Feynman prescription.
As already mentioned in the introduction, the amplitude has several physically interesting limits, that we discuss presently, as summarized in figure 1.

Soft limit
When |s|, |t| ≪ m 2 (keeping s/t fixed), the massive W bosons can be integrated out, leading to a local effective action. At tree-level, the massive W bosons do not appear when scattering the light SU(4) particles, so that the scattering amplitude is the same as in the unbroken theory. On the other hand, at loop level (and in the large N c limit), the light particles do not interact directly among themselves, but through a loop of massive W bosons. We have In this formula, the g 2 /(6m 4 ) term is one-loop exact, as predicted from non-renormalization theorems (see ref. [31] and references therein).

Forward limit and total cross-section
In the forward limit t = 0, the optical theorem relates the imaginary part of the scattering amplitude of massless scalars YȲ −→ YȲ to the total cross-section of Y,Ȳ producing a JHEP04(2018)047 pair of massive W bosons, plus other particles. We have [32] where E cm = √ s is the center of mass energy and p cm = √ s/2 is the center of mass momentum of one particle. We have In the Euclidean region −s > 0 we find Analytically continuing to s > 4m 2 (taking into account the Feynman i0 prescription), and taking the imaginary part of eq. (2.10), and using formula (2.8), we arrive at In appendix A, we compute this cross-section to the three-loop order, and observe a simple pattern, which we argue allows us to propose an exact formula of its high energy limit: (2.12) Here B(g 2 ) is the Bremsstrahlung function, given below in eq. (2.21) as an exact function of the coupling. It is striking, in particular, that the total cross-section for massless scalars or photons remains finite as s → ∞. This is likely a consequence of working at the leading order in the large N c limit (which neglects, in particular, the interactions within the unbroken SU(4) subgroup). It is instructive to recall Pomeranchuk's theorem [33], which states that if the crosssection grows with energy, the cross-section for a particle and its antiparticle will be asymptotically equal. The hypothesis of the theorem is not satisfied: here the cross-section does not grow with energy. Interestingly, the conclusion is also maximally violated: the amplitude for the antiparticle process vanishes at this order:

Threshold expansion
Let us consider the amplitude close to the threshold s = 4m 2 for producing a pair of W bosons. We expect the perturbative series to break down in the regime when the velocity β u ∼ g 2 , for the following reason. The produced W bosons interact by exchanging massless gauge fields and scalars (and fermions), from the unbroken part of the gauge group, which lead to an attractive 1/r potential. This causes the W bosons to form non-relativistic Hydrogen-like bound states, which are exactly stable in the large N c limit. Their binding energies are of order ∆E ∼ mg 4 at weak coupling. While one cannot see these bound states in our fixed-order calculation, one should still expected to see the perturbative series JHEP04(2018)047 diverge when the kinetic energy becomes of this order. Recalling that β u = 1 − 4m 2 /s, this indeed translates to β u ∼ g 2 .
Physically, the leading terms should originate from a nonrelativistic hydrogen-like system with the Hamiltonian in the center-of-mass frame 2 The contribution of this system can in fact be computed analytically as a function of g 2 /β u (see ref. [34], eq. (4.55)): (2.14) To all orders in g 2 this predicts the most singular term as the velocity β u → 0. This resummation accounts for certain ladder graphs; these are the same graphs which govern the Regge limit. There is in fact a very close connection between these two limits, as discussed in ref. [12]. Higher order corrections to eq. (2.14) should be interpreted as relativistic and many-body corrections to the Coulomb Hamiltonian.

High energy limit
We can take s, t to be much larger than the mass, m 2 /s → 0, m 2 /t → 0, with s/t fixed.
In this case, the mass serves as an infrared and collinear regulator, leading to double logarithms of the small mass. Expanding eq. (2.4) in this limit, we obtain (2.15) It was argued [11], based on anomalous dual conformal Ward identities originally derived for Wilson loops [10], that the four-point amplitude should have the following exact form, where γ(a) is the light-like cusp anomalous dimension [35,36],G 0 is a collinear anomalous dimension, andc a coupling-dependent constant. Eq. (2.16) can be viewed as a massregulated version of the BDS ansatz [9], which was originally formulated within dimensional regularization. The small mass limit and eq. (2.16) were studied previously using Mellin-Barnes techniques in refs. [21,22]. In the preceding sections, we derived analytic formulas for M up to three loops. As a check, we reproduced eq. (2.16) to that order by taking the small mass JHEP04(2018)047 limit of our formulas. For reference the coefficients are γ(g 2 ) = 8g 2 − 16ζ 2 g 4 + 176ζ 4 g 6 ; G 0 = −4ζ 3 g 4 + g 6 (36ζ 5 − 8ζ 2 ζ 3 );c(g 2 ) = 3g 4 ζ 4 − g 6 (50ζ 6 + 16ζ 2 3 ). The fact that the logarithm of the amplitude does not grow faster than log(s) is consistent with behaviour expected in the Regge limit s ≫ m 2 , that will be discussed presently. The fact that eq. (2.16) contains only a single logarithm (and no further s dependence) means that M is Regge exact (in the small mass limit). 3

Regge limit
The leading term in the Regge limit s ≫ m 2 , t, up to power corrections, has been discussed in refs. [21,22]. It is given by a single power law, The leading terms are given by We note that in planar N = 4 super Yang-Mills, the Regge trajectory is equal to the angle-dependent cusp anomalous dimension [21], The angle-dependent cusp anomalous dimension is an extremely interesting quantity in its own right. In QCD it is known up to three loops [39], while in N = 4 sYM theory it was computed up to three loops in [26] and in the planar limit up to four loops in [40]. Furthermore it is controlled by integrability [7,8] in N = 4 sYM theory. In particular the method of ref. [46] can be used to evaluate j 0 numerically at any value of the coupling. We refer the interested reader to [26,39] for a discussion of its various properties. Here we wish to point out that its small angle limit is known exactly [41], Here λ = g 2 YM N c , and I 1 and I 2 are Bessel functions. The φ 4 -term in the small angle expansion of Γ cusp (φ) is also known analytically [46].
Computing the three-loop Regge limit (2.17) of the amplitude, as described in appendix C, we computed j 0 (t), and hence Γ cusp (φ), to the three-loop order. In this way, we reproduced the result of ref. [26]. For the explicit expressions we refer to that reference.

Discussion
In summary, we studied a scattering amplitude of four complex scalars in N = 4 sYM theory with a spontaneously broken gauge group, which is known as an analytic function of the variables u = 4m 2 /(−s) and v = 4m 2 /(−t) up to three loops [1]. We studied several kinematic limits. To obtain these limits we used techniques for differential equation and calculated the amplitude in an (asymptotic) expansion, in principle to any order in the expansion parameter. For convenience of the reader, we provide our results as computer-readable supplementary material. Interestingly, we find that the leading terms of most of the above limits are in principle known to all loop orders, or are controlled by an integrable model.
In the case of massless amplitudes, a systematic expansion around the collinear limit could be found and described via integrability [42,43]. The above observations nurture the hope that something similar can be done for massive amplitudes, at least in one of the above limits. Focusing on the first subleading term in the Regge limit, we will find a remarkably simple structure. In the remainder of this paper, we discuss these findings in detail.

Regge expansion using dual conformal partial waves
The Regge limit is a likely candidate around which one can hope to build a systematic expansion. When discussing power-suppressed terms, however, the choice of variable used to express the leading term (2.17) becomes important. In this section we will use symmetries to derive a good parametrization of the limit, and we will find that the first power-correction is then controlled by a single independent power.

Dual conformal partial wave analysis
A simple improvement over expanding in inverse powers of 1/s in the Regge limit at fixed t is to use instead the partial wave expansion, where each power gets upgraded to a Regge pole. Each Regge pole contributes proportional to a Legendre function, producing an asymptotic expansion of the amplitude in the form (see eq. (A.12) of [44] or (2.9.6) of [45]): where j n (t) denote the Regge trajectories and Q −j−1 are associated Legendre functions.
In the Regge limit these tend to 4 This highlights how each Regge pole resums an infinite tower of inverse powers of s. The Legendre functions which control each tower originate physically from the SO(2,1) Lorentz subgroup which preserves the spacelike momentum exchanged in the t-channel. For this reason, the expansion (3.1) generally contains fewer independent coefficients than 4 The right-hand side shows the expansion of: − tan(πj)Γ(j+1)

JHEP04(2018)047
the straightforward 1/s expansion. This is analogous to how conformal blocks are used to resum descendant operators in conformal field theories. The N = 4 sYM amplitude that we consider benefits from a larger SO(4) or SO(3,1) dual conformal symmetry. It was identified in [12] as a relativistic version of the Laplace-Runge-Lenz symmetry of the hydrogen-like bound states which appear in intermediate states of the amplitude. This symmetry implies further relations among the Regge trajectories.
Using the embedding formalism to work out the action of the dual conformal symmetry on the kinematical invariants, and the detailed form of SO(4) Legendre functions, a symmetry-improved expansion around the Regge limit is derived in appendix B: where the variable, which vanishes in the Regge limit s → ∞, is (recall eq. (2.6) for our notations) Again each term behaves in the Regge limit like a power of s: Y −j−1 ∝ s j+1 , however an intricate tower of subleading powers, similar to but distinct from eq. (3.2), is associated with each SO(4) Regge trajectory j n (t). The angle (3.4) is distinct from the usual scattering angle between the two external massless photons (see eq. (3.1)), due to the −s/(2m 2 ) term. This makes the angle real in the physical region of t-channel scattering when 0 < t < 4m 2 and −t < s < 0. It is hyperbolic above the massive continuum t > 4m 2 , or whenever |s| is sufficiently large, as is assumed in the Regge limit formula (3.3).
As an application of this formula, we have subtracted from the three-loop amplitude described in the preceding section its leading behavior: r 0 (t)Y −j 0 (t)−1 , which is known to exponentiate to a single power law (with exponent previously given in refs. [26,46]). Amazingly, looking at the remainder, we find that the first subleading power also exponentiates! That is, we find to three loops that the logarithm of the remainder is linear in log Y . We can thus write where j 0 ≈ − 1 and j 1 ≈ − 2. This is rather remarkable: in principle the first subleading term could have been a sum of multiple power laws. In fact, had we used the straightforward 1/s expansion, or the usual SO(3) Regge pole expansion, we would have (artificially) found at subleading order two distinct power laws, with one exponent mysteriously equal to the leading exponent minus one, thus hinting at the presence of a hidden symmetry.
Since we have only one single subleading power, taking the limit of the three-loop amplitude of ref. [1] using the method detailed in appendix C, allows us to extract its exponent to three loops. Defining

JHEP04(2018)047
our result for the Regge trajectory is: where H are harmonic polylogarithms [47,48] with argument 1 − x 2 . The residue starts as We also looked at the next powers in the expansion. We find that there is never any double logarithm, that is, never more than one power of logarithm per loop order. To leading logarithm order, the sub-sub-leading term in the amplitude is given as One can show that this not the exponential of a single power. For convenience, the full three loop expansions of j 0 , j 1 , r 0 , r 1 and the sub-subleading amplitude are recorded in the supplementary material attached to the arXiv submission of this paper.

From the Regge limit to Wilson lines with a cusp
While we currently lack a systematic effective field theory framework to characterize the subleading powers (3.5) in the Regge limit, to make progress we will use a map to an equivalent problem involving a cusped Wilson line, following [12,21,26]. The idea is to generalize the Higgs symmetry breaking pattern by further breaking SU(4) down to U(1) 4 , thus allowing a distinct mass m i for the W bosons. The amplitude still depends only on six-dimensional dot products of the vectors (B.2), which give rise to two independent cross-ratios, now equal to [49,50]: (3.9) These reduce, in the equal-mass case, to our previous definitions, e.g u = 4m 2 −s , see figure 2 (a). That the amplitude M depends only on these two cross-ratios has an implication which is familiar in the non-relativistic limit: the cross-ratios then depend only on the kinetic energy divided by the reduced mass. Most important for us, will be the fact that the Regge limit s → ∞ of this amplitude, is equivalent to the massless limit m 3 → 0. See  Figure 2. The amplitude on the left, with four massive W bosons (thick lines) running outside the loop, is equivalent, through eq. (3.9), to an amplitude with unequal-mass bosons. In a limit equivalent to the Regge limit, one of the masses go to zero, revealing infrared divergences within the associated cusp.
Logarithms in such a massless limit are naturally associated with soft quanta coupled to Wilson lines with a cusp geometry, see figure 2 (c). This was used in [26] to obtain the 3loop cusp anomalous dimension from the leading power in the Regge limit of the amplitude. It was used in the other direction in [12], to obtain the bound state spectrum of hydrogenlike states from the latter (and further analyzed at strong coupling [51]). Using the same map for the first subleading power, the subleading Regge trajectory j 1 , see eq. (3.5) can thus be used to predict that a Wilson line with a cusp has an excitation with scaling dimension: where j 1 is given to three loops by eq. (3.7). The variable x used in that expression can be expressed in terms of the cusp angle φ between p 3 and p 4 shown in figure 3 as In summary, the fact that the first power correction in the Regge limit exponentiates is remarkable and strongly suggests the existence of a systematic expansion with a nice structure. In the next section, we make a proposal for how to determine j 1 (t), or equivalently the scaling dimension Γ cusp,Φ (φ), in terms of an effective field theory calculation.

Renormalization of Wilson lines with operator insertions
In this section we show perturbatively up to two loops that the exponent j 1 in the power suppressed term in the Regge limit can be identified as the scaling dimension of a cusped Wilson loop with a scalar inserted at the cusp point.

Cusped Wilson line
We recall that we found for the exponent j 0 of the leading term in the Regge limit the relation  where Φ is one of the six scalars. The fields are taken to be in the adjoint representation of the gauge group. 5 In the following, we will work in the large N c limit. The contour consists of two straight line segments which form a cusp at the origin. Note that the directions are chosen such that v 1 is incoming and v 2 is outgoing; see figure 3. The cusp angle is defined as The Wilson loop operator we are interested in is then simply given by  [53][54][55][56][57][58][59]. In a conformal field theory, such as N = 4 sYM, the renormalization factor Z cusp has the following form in D = 4−2ǫ dimensions, see e.g. [60], cusp is the (angle-dependent) cusp anomalous dimension. To motivate which subleading power operators to consider, let us discuss schematically what might be anticipated from a systematic analysis of the massless limit m 3 → 0 of the preceding subsection, in the framework of heavy quark effective theory (HQET). One would start by integrating out the heavy (internal) W bosons, which should reduce the 2 → 2 amplitude at leading power to a matrix element of local operators which creates a pair of heavy (scalar) quarks: Taking the matrix element in a state with two heavy particles, the propagator for the HQET fields q simply produce Maldacena-Wilson lines, thus recovering W cusp . Schematically, two kinds of power corrections are then expected: either from higher-dimension operators in the HQET Lagrangian, or from higher-dimensional corrections to the local operator. The former adds operator insertions along the Wilson lines, whereas the latter adds fields or JHEP04(2018)047 derivatives strictly at the cusps. For the present analysis, we shall assume that the leading power corrections come entirely from cusp insertions, and ignore the corrections to the Wilson lines.

Operator mixing and renormalization
For the first power correction, we are led to consider local HQET operators with mass dimension one higher than (4.5), for example q † Φq. As we ignore corrections to the HQET Lagrangian, the heavy quark fields become again simply Maldacena-Wilson lines, and we can treat this as a scalar insertion in W cusp . Since the small mass m 3 is controlled by the Higgs mechanism, it should be better viewed as a property of the state rather than of the operator, and in the free theory the scalar Φ simply becomes its expectation value m 3 . The considered operator can in principle mix with any other which has the same mass dimension and Lorentz indices. The only gauge invariant operator built from N = 4 sYM fields fulfilling these criteria involve derivatives within the plane of the cusp (and therefore total derivatives). We are thus led to consider the two sYM operators: The relative sign in (4.6) comes from the symmetry under v 1 ↔ −v 2 .
In the case where the scalars which couple to the Wilson line are orthogonal to the scalars inserted at the cusp, the anomalous dimension is known from integrability [46], see also [61] for related work using localization techniques. Here, however, these scalars are the same. This setup was considered by Alday and Maldacena in [62], where they computed the anomalous dimension at one loop for the straight line case. We extend the calculation to two loops with a dependence on the cusp angle.
In general both operators can mix at loop level. In order to resolve this mixing, we will consider suitable correlation functions of these operators. Since 0| W cusp,Φ |0 has no treelevel contribution, we find it more convenient to consider the correlator with an additional scalar Φ(p 3 ). We take the latter to be on-shell, p 2 3 = 0, so that the correlators are gauge independent.
At tree-level, we find These correlators depend on the cusp angle φ, as well as on the two invariants s 1 = −2v 1 ·p 3 and s 2 = +2v 2 ·p 3 . In the following, the different momentum-dependence of the correlators, together with the fact that the (ultraviolet) renormalization matrix must be independent of the external momentum p 3 , will allow us to resolve the operator mixing. Moreover, since ∂W cusp is a derivative of the lower dimensional operator W cusp , it renormalizes multiplicatively with the same renormalization factor, and we thus expect JHEP04(2018)047 the mixing matrix to be triangular. Taking also into account that operator mixing only appears at the loop level, we have the following structure of the renormalization matrix for In the following, we show up to two loops that Z is diagonal, i.e. that W cusp,Φ renormalizes multiplicatively up to that loop order. This accidental vanishing of Z mix might be related to the enhanced (dual conformal) symmetry of our setup, which we are not exploiting in the present calculation.
The correlators not only have UV divergences coming from the cusp and the operator insertion, but also soft and collinear divergences from the on-shell scalars. The latter can be renormalized with a common IR Z-factor where µ is the renormalization scale, and γ (L) is the coefficient of (g 2 ) L in the cusp anomalous dimension defined below (2.16); γ HgH , discussed below, is interpreted physically as the SCET collinear anomalous dimension of one massless field, plus that of two heavy fields. This Z-factor simultaneously renormalizes the IR divergences of both correlators, because the structure of the IR divergences arises only from the configuration of the external lines and not from the cusp point. Likewise UV divergences only come from the cusp and the operator insertion and not from the external scalar. The renormalization condition is then given by Given the known cusp anomalous dimension (4.4), this equation allows us to determine the IR Z-factor and then the missing pieces in Z −1 .

One-loop calculation
Let us now discuss in some detail the one-loop calculation. In figure 4 and 5 some sample Feynman diagrams are shown. We find the following result for the correlators, The IR Z-factor (4.12) is determined by the first component of (4.13), with the result The first matches the expansion of the cusp anomalous dimension γ(g 2 ) = 8g 2 − 16g 4 ζ 2 + O(g 6 ) stated earlier.
The second component of (4.13) then allows us to calculate the remaining pieces of Z −1 (or, equivalently Z). We can deduce that at the one-loop level, the lower off-diagonal element in Z −1 has to vanish. This is seen as follows. The tree level contribution of (4.14) to the second component in (4.13) has the form [g 2 Z For the general case of a renormalization matrix, the corresponding anomalous dimension matrix Γ is given by Here we used that the β-function vanishes in N = 4 sYM theory. In our case the renormalization matrix is diagonal, and therefore the matrix inversion is trivial. In this way, adding the engineering dimension of the scalar insertion, we obtain This is in agreement with our conjectured relation (3.10) (and also with the computation of the same anomalous dimension in [62]).

Extension to two loops and discussion
We proceeded to perform the calculation to two loops. Some details of this calculation, and in particular the computation of the necessary Feynman integrals, can be found in appendix D. Proceeding as in the one-loop case and subtracting the known ultraviolet divergences of W cusp , we find for the IR renormalization coefficients The cusp anomalous dimension matches the expected result, but it is also instructive to compare the value of the second, single-logarithmic infrared divergence. According to ref. [63], it should be the sum of a collinear anomalous dimension for massless particles, plus the infrared anomalous dimensions for the two massive fundamental Wilson lines (or equivalently heavy particles): γ HgH = γ g + 2γ H,fund . (4.20) The known two-loop collinear anomalous dimension for any parton in this theory is γ g = 2g 4 ζ 3 , which matches with the maximal transcendental part of the QCD result for either a gluon or an adjoint quark (see [64]). The massive case γ H,fund was not calculated, to our knowledge, however we can take the maximal transcendental part of the result for an adjoint QCD quark as given in [63]: γ H,adjoint = −4g 4 ζ 3 ≡ 2γ H,fund . These indeed sum up to the value (4.19), which we conclude is consistent with the principle of maximal transcendentality and the known QCD values.
Proceeding to the UV renormalization, we find that up to two loops, Z −1 remains diagonal, and gives the anomalous dimension where H are harmonic polylogarithms with argument 1 − x 2 . This is in perfect agreement with the conjectured relation (3.10) with the Regge trajectory given in eq. (3.7). This two-loop calculation provides non-trivial evidence in favor of our proposed relation (3.10), which identifies the subleading Regge trajectory j 1 with the scaling dimension Γ cusp,Φ of a Wilson loop operator. We note that the Regge trajectory j 1 is known (from the scattering amplitude calculation) to one more order, three-loops. It would be interesting to verify that eq. (3.10) also holds at that order, perhaps using integrability or other methods.

Conclusion and outlook
In this paper we studied a massive four-particle scattering amplitude in N = 4 super Yang-Mills. Starting from the three-loop result obtained in [1], we initiated a systematic analysis of this amplitude as a function of all kinematic invariants. While the full amplitude involves complicated multiple polylogarithms that depend on two variables s/m 2 and t/m 2 , it simplifies considerably in various physically interesting limits. In appendix C, we explain in detail how we obtained expansions from the differential equations of [1], and we provide the detailed three-loop formulas for the asymptotic expansions described in section 2 as supplementary material.

JHEP04(2018)047
With the results of the asymptotic expansions at hand, we focused on the Regge limit. Using the fact that the amplitude is governed by a hidden symmetry, dual conformal symmetry, we derived a partial wave expansion which incorporates this additional information. It was known that the leading term in the Regge limit exponentiates, with the exponent given by the anomalous dimension of a Wilson loop with a cusp. In this paper, we explored subleading, power-suppressed terms. Surprisingly, we found that, in this expansion, the first power suppressed term also exponentiates! We computed the Regge exponent j 1 to three loops, cf. eq. (3.7).
Moreover, by using the symmetry to map the Regge kinematics to a soft expansion, cf. section 3.2, we argued that the relevant operator controlling the subleading Regge limit should be a cusped Wilson loop with a scalar insertion at the cusp. We verified this proposal to two loops in perturbation theory. In order to do so, we performed a soft current calculation for massive quarks, to two loops, and found a perfect match. The details of the calculation of the Feynman integrals are presented in appendix D.
While we found that all integrals needed to compute the divergent part of the two-loop correlators required only multiple polylogarithms. It is interesting to mention that some of the finite integrals involve elliptic polylogarithms.
We briefly discuss some interesting questions for future work. Our asymptotic expansions provide a wealth of data to explore higher order terms in the power expansion. At subsubleading power in the Regge limit (1/s 2 ), we compared the expansion (3.8) with an ansatz with two power laws. We find that such an ansatz is inconsistent with angle-independent one-loop exponents, a property which would be expected in the perturbative expansion. However, it is possible to write a consistent ansatz with three exponents. Such an ansatz could be tested once higher-loop results for the scattering amplitude become available.
It would be interesting to derive a systematic expansion using heavy quark effective theory as sketched around eq. (4.5), exploiting the setup of section 3.2 where the problem is mapped to a massless limit m 3 → 0. The leading order terms correspond to the cusped Wilson loop, while we proposed that at first order one needs to consider a scalar insertion into the Wilson loop. However, in general, the Lagrangian of this effective theory contains an infinite series of irrelevant operators suppressed by the heavy mass, giving corrections to the Wilson lines, analogous to those considered in [2,3,65,66]. At higher orders in the power expansion, these may cause triangular mixing between the operator (4.5) and higherdimensional ones. This could prevent the amplitude from being a sum of pure power laws. It would thus be very interesting to elucidate the structure at higher powers, and ultimately to translate these findings to the usual null Wilson lines approach to the Regge limit.
Once the operators are identified, a separate question consists in computing their anomalous dimensions. The cusp anomalous Γ cusp is known to be governed by an integrable system [7,8], which was simplified in [67]. It would thus be interesting to see if the first subleading trajectory, given to three loops in eq. (3.7), can be reproduced quantitatively by extending the methods used in [46]. We also wish to point out that though the AdS/CFT correspondence, it may be possible to study the anomalous dimensions at strong coupling [68].

A The total cross-section to three loops
Here we discuss the total cross-section σ tot = σ YȲ →X , see eq. (2.8), up to three-loop order, using the results of [1]. We will observe an interesting property about its high-energy behavior that motivates the all-loop prediction (2.12).
In [1] the integrals contributing to the three-loop amplitude were computed. They fulfill a differential equation with trivial boundary conditions at s = t = 0. To study the forward limit t = 0 we thus only need to integrate with respect to s, or equivalently u = 4m 2 −s . The differential equations have logarithmic-type singularities at u = 0, −1, and a squareroot type singularity at u = −1. Upon switching to the variable x = βu−1 βu+1 , the alphabet becomes {log(x), log(1+x), log(1−x)}. From this it follows that the limit of M can be written in terms of harmonic polylogarithms [47,48] with argument x. Our result takes the form The two loop result originates solely from the horizontal ladder, s 2 tG 1,1,1,0,1,0,1,1,1 , because the other integral is explicitly proportional to st 2 . For the same reason, at three-loops only two integrals from [1] contribute to the forward limit.
The first term originates from the ladder while the last line originates from the tennis court diagram. Here H are harmonic polylogarithms of argument x, which we omitted for brevity. Note that all formulas above are symmetric under x → 1/x, as may be verified by using identities between the harmonic polylogarithms for different arguments [47,48].
The physical region above the threshold s > 4m 2 , where the inelastic cross-section (2.8) is nonzero, correspond to −1 < x < 0. This region can be reached by analytic continuation. In total we can identify three regions relevant for the analytic continuation: the Euclidean region, below the threshold and above the threshold. In the Euclidean region and below the threshold the amplitude is real. Above the threshold, however, the amplitude has a branch cut and following Feynman's i0 prescription we have to evaluate the harmonic polylogarithms slightly above the real axis. The analytic continuation is illustrated in figure 6. By calculating the imaginary part of the amplitude we find for the total cross-section where Here we have chosen a form that is manifestly real-valued for −1 < x < 0. The cross-section (A.4) can be seen to approach a constant in the high-energy limit x → 0: Remarkably, this agrees precisely with the perturbative expansion of the Bremsstrahlung function (2.21)! This is not a coincidence. To see this, we use the leading Regge behavior

JHEP04(2018)047
of the amplitude at large s and fixed t wherer 0 and j 0 + 1 are given in (2.18). Because of the mass gap of the W bosons, loop corrections to the amplitude A ∝ M s/t must be real and analytic around t = 0, which implies that the coefficient of 1/t is tree-level exact. Thus the loop corrections to these parameters must vanish at the origin:r 0 (0) = 1 and j 0 (0) = −1. Furthermore the coefficientr 0 is real, so the imaginary part originates from the trajectory. Thus Equation (A.10) is a bit unusual since the cross-section involves the slope of the Regge trajectory at t = 0 rather than the intercept, as is more usual. This happens here because the intercept precisely vanishes. Using the relation (2.21) expressing the slope at t = 0 in terms of the (exactly known) Bremsstrahlung function gives the prediction (2.12) for the total cross-section, generalizing (A.8) to all orders.

B Regge expansion using dual conformal partial waves
As discussed in section 3, the enhanced symmetry of the amplitude we look at makes it possible to efficiently organize the Regge limit. In this section we derive this improved expansion which exploits dual conformal symmetry. It will be helpful to use the embedding formalism, which realizes Minkowski space as the null cone in R 4,2 . The region momenta associated to each (planar) loop are represented as a (projective) 6-vector where the first four components (before the vertical line) are spacelike and the last two are timelike. Here µ is an arbitrary scale and x 2 = x 2 − (x 0 ) 2 . In the presence of internal masses we also need to consider timelike dual coordinates for the external regions (see figure 7), satisfying Y 2 i = −m 2 : These definitions ensure that six-dimensional dot products give massless and massive momentum space propagators: This notation is helpful because the dual conformal symmetry SO(4,2) acts linearly as rotations of these 6-vectors. The external momenta of our planar four-particle scattering problem are encoded in the differences between the above points Y A i as: definiteness let us begin by assuming kinematics with a timelike t-channel, with 0 < t < 4m 2 where t = −(p 2 + p 3 ) 2 . We can use Lorentz invariance to go to the rest frame of p 2 + p 3 and use translation invariance in y-space to set y 2 = y 4 = 0. Furthermore the energies of p 1 , p 2 must be equal and opposite, which allows to set y 0 1 = y 0 3 = 0. In this frame the dual coordinates reduce to where α = √ 4m 2 − t and in addition we have chosen µ = α 2 in order to set to zero the last spacelike component of Y 2 and Y 4 .
The SO(4) symmetry is apparent in this frame: these are simply the rotations of the first four components, which preserve the two dual coordinates Y 2 and Y 4 bounding the t-channel. This contains the usual SO(3) rotations of a pair of particles in its rest frame, with three additional generators which are related, in the nonrelativisitc limit, to the flow generated by the Laplace-Runge-Lenz vector of the hydrogen atom [12].
From eq. (B.4) we see that the dependence on the Mandelstam variable s = −(p 1 +p 2 ) 2 is encoded in the SO(4)-invariant angle between the first four components of Y 1 and those of Y 3 . From a short computation: where we have simplified using that p 2 1 = p 2 2 = t 4 and − p 1 · p 2 = t 4 + s 2 . As noted in the text, this angle differs from the usual scattering angle between the two external massless photons (see eq. (3.1)) by the −s/(2m 2 ) term. It is real, for example, in the t-channel region where 0 < t < 4m 2 and −t < s < 0.
We now derive the corresponding partial wave expansion, starting from the case where the angle is real and then analytically continuing. The idea is to express the dependence on s in terms of SO(4) spherical harmonics; for each SO(4) spin j these sum up to a Chebyshev polynomial (the SO(4) analog of the Legendre polynomials):

JHEP04(2018)047
Using the relation similar to (2.2) between the YȲ →Ȳ Y amplitude and the stripped matrix element M , and absorbing s-independent factors into the coefficients c j (t), we can rewrite this as an expansion for M : We note that, as a mathematical statement about bases of functions, one could equally well apply this decomposition to M itself (or to M times any function of the cross-ratios u, v). The amplitude A is singled out physically since its t-channel cuts have a Hilbert space interpretation in terms of intermediate states acted upon by SO(4). 6 We now have an expansion valid for real angles. To reach the Regge limit Y → 0 where the angle is imaginary, we follow the standard procedure and rewrite the sum as an integral using the Watson-Sommerfeld trick [44,45]: where we have assumed that 0 < θ < π and the phases have been chosen such that the integrand vanishes at large imaginary j (assuming that c j (t) is bounded). The idea is that, deforming the contour to the right and picking up the residues of sin(πj), this reproduces the sum (B.8). But taking the Regge limit Y → 0, a different contour deformation becomes appropriate. The contour can still be closed to the right in the first term, but now to the left in the second term. Two types of singularities arise: poles from the inverse sine factor, which add up to: These however neatly cancel out, because for integer spin j the coefficients c j are odd under j → −j −2. Such a cancellation of kinematic poles can be proved from the Froissard-Gribov inversion formula and occurs generally for any SO(D) expansion; a detailed discussion in the SO(3) case can be found in refs. [44,45]. All that remain are the Regge poles of c k (t) from the second term, which add up to the asymptotic expansion: where we have defined the residues r n (t) = πe iπjn(t) sin(πjn(t)) Res j=jn(t) c j (t). This formula is used in the main text to efficiently organize the Regge expansion Y ∼ 1 s → 0. JHEP04(2018)047 high energy limit s, t → ∞ Figure 8. Different limits we consider in the u − v plane. To derive expansions, first the boundary value for each limit is obtained. Initially known in the soft limit, the boundary value is transported along the edge of the diagram.

C Method for expanding the three-loop amplitude of [1]
Here we explain how to express the amplitude in various limits, which in general can contain logarithmic divergences. In principle, we could use the analytic expressions for the master integrals derived in ref. [1], and expand them using properties of the iterated integrals they were expressed in. We find it more convenient to obtain such expansions directly from differential equation for the master integrals that were derived in ref. [1]. In order to do so, we use a well-known procedure for solving differential equations in a limit, following closely the textbook [69]. Let x be parameter that parametrizes the expansion around x = 0, and let f be the vector of master integrals. As we will see, the solution for f takes the general form P (x)x A 0 f 0 , where P (x) is a (matrix) polynomial in x; the matrix exponential x A 0 contains possible logarithmic divergences, and f 0 is the finite boundary value at x = 0. Given possible powers of logarithms log(x), one may also call f 0 the 'regularized' boundary value.
A technical point is related to obtaining such boundary values for all expansions that we are interested in. The boundary value considered in ref. [1] is taken at s, t → 0, see figure 8. In order to obtain appropriate boundary values for other expansions, we first transport this value to other regions, along appropriate paths. By 'transporting' we mean solving the differential equation along a given path. In principle one could choose any convenient path. However, some choices are preferable over others. In particular, one can often find paths for which the one-parameter solution is expressible in terms of a relatively simple class of functions, the harmonic polylogarithms. This is the case for the paths shown in figure 8.
As we will discuss in more detail in the following, special care is required when singular boundaries are approached (corresponding to singularities of the differential equation). When several of such boundaries intersect, it is important to clarify how the singular boundary is approached. In mathematical language, one can perform a 'blowup' that resolves singular intersections of boundaries.

JHEP04(2018)047
As a non-trivial verification of our analytic continuation procedure, we verified that, upon returning to the original point s, t → 0 after going around the whole square in the positive quadrant shown in figure 8, we recover the correct boundary value.

C.1 Solving the differential equation in an expansion
In this section we follow [69] closely. Given a square n-th order matrixĀ(x), which is holomorphic on a connected open set R ⊂ C, the differential equation f ′ (x) =Ā(x) f (x) has a unique solution on R, provided a boundary condition f (a) = f BV , a ∈ R. Furthermore this solution is holomorphic on R. We are interested in the more special case where the matrixĀ(x) has a regular singular point x p / ∈ R. Without loss of generality we choose x p = 0. Then the differential equation can be rewritten as As a first step in solving (C.1) we perform a transformation f (x) = P (x) g(x) with a non-singular holomorphic matrix P (x) to get The new matrix B(x) is determined by P (x) and A(x). Our aim is now to find a P (x) such that B(x) becomes as simple as possible, in order to solve (C.2). As we will see, in practice, we can choose B(x) and then calculate P (x) using Inserting the respective power series for A(x) = k∈N 0 A k x k as well as for B(x) and P (x) in the differential equation above, we obtain, after equating the coefficients, the following recursion relation At this point a subtleness arises: we are of course interested in an unique solution of the problem, but the matrix equation AX − XB = 0 for given square matrices A and B can in principle has a non-trivial solution for the matrix X. One can show that the equation AX − XB = 0 has such a non-trivial solution X = 0 if and only if A and B have at least one common eigenvalue. Equipped with this knowledge, we now choose a certain matrix B(x). If A in (C.1) is a constant matrix, then we do not need to simplify the problem any further. Therefore we choose as our starting point B 0 = A 0 and P 0 = 1.
From the previous argument we know that if no pair of eigenvalues of the matrix A 0 differs by a positive integer, the P k are determined by (C.4). In our case the matrix appearing in the differential equation for the master integrals is a lower triangular matrix with vanishing diagonal elements, hence all eigenvalues are zero and the former condition JHEP04(2018)047 is trivially fulfilled. We choose B k = 0 for k > 0 in order to obtain a simple differential equation after the transformation. With this choice B(x) = A 0 the solution of (C.2) is given by x A 0 . Returning to the original problem, we find the asymptotic expansion of the solution of (C. 1) f where P (x) is calculated recursively from (C.4) using P 0 = 1, B 0 = A 0 and B k = 0 for k > 0. We wish to make the following comments.
• As was already mentioned earlier, the solution may have logarithmic divergences in the limit x → 0. If present, these are described by the matrix exponential exp[A 0 log(x)]. We call f 0 the boundary value at x = 0, even in such singular cases.
• One can interpret the matrix as the fundamental system of solutions of the matrix differential equation In the following subsections, we describe in more detail the procedure of obtaining the boundary values for the different expansions, and on the choice of variables for the latter.

C.2 Soft expansion
The soft or low energy limit describes the region where |s|, |t| ≪ 4m 2 . We perform the calculation in the Euclidean region s, t < 0, but the result is valid in the entire region |s|, |t| ≪ 4m 2 , since the result is simply a polynominal in s and t.
In order to derive the soft expansion, we can use our starting boundary value at m → ∞ with g start = (1, 0, . . . , 0). For solving the differential equation we introduce the transformation with the ratio R = u/v = t/s fixed. The differential equation is then solved for small x, as explained above. Applying this to the amplitude, we find, up to three loops, 24 . (C.7) As mentioned earlier, the 1/m 4 term is one-loop exact. In the perturbative expansion this feature is evident already pre-integration: all higher-loop integrals appearing in the perturbative expansion are explicitly proportional to at least s 2 t or st 2 , e.g. see figure 7 of [70] for the five loop integrand. It is noteworthy that all coefficients are rational multiples of g 2 = g 2 Y M N/(16π 2 ); no transcendental numbers such as ζ values appear. Technically this can be traced to the fact that in ref. [1] a uniform weight basis could be found in which all but one integral vanish in the low-energy limit. It would be interesting to see if this remains the case at higher loop orders.

C.3 Regge expansion
Switching from the kinematic invariants s and t to the variables u and v, the Regge limit is described by v ≫ u. In our calculation we consider the limit u → 0 with v > 0. To transport the boundary value from our starting point at (u, v) = (∞, ∞) (this is the limit m → ∞ in the Euclidean region) to our end point, we split the path in two straight line segments γ 1 and γ 2 . The first path γ 1 = {(u, v) = (−t, ∞)|t ∈ (−∞, 0]} is parallel to the u-axis and ends on the v-axis, while the second path γ 2 = {(u, v) = (0, −t)|t ∈ (−∞, −ṽ]} is on the v-axis and ends at someṽ > 0.
For the first segment γ 1 we take the limit v → ∞ and substitute the variable u: This leads to a differential equation for the master integrals g on the path γ 1 with the Since A γ 1 is a lower triangular matrix we obtain the solution recursively. It can be expressed in terms of harmonic polylogarithms and constants, which are determined by the boundary value at x u = 1. For the path γ 2 we need the boundary value at x u = 0. Unlike at x u = 1 the master integrals exhibit logarithmic divergences at x u = 0. To extract the (regularized) boundary value, we use the asymptotic expansion of the differential equation.
The calculation for the second path γ 2 is identical to the previous one. In the limit u → 0 and with the variable transformation (C.8) for the variable v the differential equation can be solved on the path γ 2 in terms of harmonic polylogarithms. With the previous calculated boundary value at x u = 0 the integration constants are fixed. With the solution g γ 2 (x v ) on the second path as our boundary value we finally solve the differential equation in an asymptotic expansion with arbitrary v > 0 in x u near x u = 0 to obtain the master integrals in the Regge limit.

C.4 High energy expansion
In the high energy limit we have |s|, |t| ≫ 4m 2 . We will work in the Euclidean region. Obtaining the appropriate boundary value at (u, v) = (0, 0) requires some care, as we discuss presently. The subtlety originates from the singularity structure of the differential equation in the limit u, v → 0. This can be immediately understood by inspecting the alphabet of the differential equation [1], which contains the letters {log(u), log(v), log(u+v), log(u 2 − 4v), log(v 2 − 4v)}. It is sufficient to study these "simple" letters, because the other letters do not add more singularities in the vicinity of (u, v) = (0, 0). The corresponding singular lines are shown in figure 9(a). The fact that the latter intersect at (u, v) = (0, 0) implies that one has to specify how exactly this point is approached. The problem of a potential ambiguity can be avoided by switching to appropriate variables that resolve the way the singularity is approached. The variable transformations can also be understood as choosing more sophisticated paths near the origin to connect the boundary values, which we obtain by approaching the origin (u, v) = (0, 0) on the u-axis or v-axis. These two boundary values are denoted by g u BV and g v BV . In figure 9 these paths are shown. The first transformation or path resolves the ambiguity of the first three considered letters This transformation is sufficient at one-and two-loops, but not at three-loops, where the new letters {log(u 2 − 4v), log(v 2 − 4v)} first appear. Therefore we introduce a further transformation. The corresponding path γ ǫ in the δ-t-plane is shown in figure 9(b) and it is divided into five straight sections γ i ǫ , i = 1, 2, . . . , 5. The crucial point is now that after the transformation the ambiguities are resolved and we can take the limit ǫ → 0 on the separate sections. Then we solve the differential equation in this limit. For example, the first segment can be parameterized by In the limit ǫ → 0 the alphabet becomes {log(τ ), log(1 − τ ), log(1 + τ )} and therefore the solution of the differential equation is given in terms of harmonic polylogarithms. However a new subtleness arises, which we so far have not encountered. The solution exhibit logarithmic divergences at τ = 0 and τ = 1. This means we have to fix the integration constants using the asymptotic expansion at τ = 0. The boundary value at τ = 0 in the limit ǫ → 0 is g u BV from before. The differential equation on the other sections can also be solved in terms of harmonic polylogarithms in the same way. After extracting the boundary value at the end point of the last section γ 5 ǫ we precisely get g v BV , which is a strong crosscheck for our calculations. Additionally the solution on the third section γ 3 ǫ ǫ→0 = (δ, t) = (0, τ ), τ ∈ [0, 1] is our desired boundary value for the high energy expansion. The boundary value depends only on the ratio R = u/v via τ = t = 1/(1 + R).

JHEP04(2018)047
With the known identities between harmonic polylogarithms the boundary value can be rewritten such that only harmonic polylogarithms with argument R appear. Now we are finally in a position to calculate the high energy expansion. For this we transform from (u, v) to (ǫ, R) using and solve the differential equation in an asymptotic expansion in ǫ. The final result is then rewritten in ρ = ǫ/(1 − ǫ) 2 = (u + v)/4.

C.5 Threshold expansion
We now consider the threshold expansion of the amplitude. For this we solve the differential equation in an asymptotic expansion in β u = √ 1 + u around β u = 0 in the physical region −s < t < 0. As in the previous limits the transportation of boundary value is the most complicated part of the calculation. However we can use the results from the forward limit, where the master integrals were calculated in the limit v → ∞ for arbitrary 4m 2 < s. As the next step we extract the boundary value at the threshold β u = 0. So far we only have considered the case v → ∞ or equivalently t → 0 − , but we are interested in an expansion in β u for arbitrary −s < t < 0 or equivalently 1 < v < ∞. Therefore we solve the differential equation at the threshold in the new variable in terms of iterated integrals. This variable transformation is chosen such that the alphabet does not contain any square roots; it is polynomial, with the highest degree being four. For brevity we do not write it down here. The integration constants of the solution of the differential equation are fixed by the previously extracted boundary value. Note that the solution does not have any divergences at y = 0 or equivalently v = ∞. Finally we us this solution as the new boundary value for the asymptotic expansion of the master integrals in β u around β u = 0. We are only interested in the imaginary part of the amplitude, which simplifies the result significantly. The imaginary part of the amplitude depends only polynomially on t, whereas the real part contains iterated integrals.

D Soft current computation
For the one-and two-loop calculation we used the Feynman diagrammatic approach and evaluated the correlators in momentum space. The diagrams where generated with QGRAF [71] and the output further processed with a custom Mathematica code which expresses the result in terms of a number of scalar loop integrals. To treat the Majorana fermions we used the techniques described in [72]. The integral reduction to a set of master integrals was done with FIRE5 [73][74][75] in combination with LiteRed [76,77]. As a non trivial cross-check the calculation was done in covariant gauge. The gauge parameter drops out after the reduction to master integrals. In addition our Mathematica code was tested JHEP04(2018)047 Figure 10. One-loop soft current integral topology. by reproducing the two-loop cusp anomalous dimension [26,39,40] and the two-loop jet function in soft-collinear effective theory [78].
In this appendix we discuss the calculation of the one-and two-loop master integrals with the differential equation method [6,[79][80][81][82][83][84]. We explain the one-loop case in detail. As a cross-check we compared the analytic results for the master integrals on several kinematic points in the Euclidean region with the numerical results obtained from FIESTA4 [85]. We found perfect agreement within the error bars.

D.1 One-loop master integrals
Let us briefly recall the kinematics of our problem. With the two directions of the Wilson lines v 1 and v 2 , and the on-shell momentum p 3 of the external scalar, we can build three invariants where we used v 2 1 = v 2 2 = 1 and p 2 3 = 0. In the following, the variable x = e iφ turns out to be most useful. At this point it is also worthwhile to introduce the Gram determinant G(q 1 , q 2 , q 3 ) = det(q i · q j ) formed by the three vectors. It is given by If all the vectors lie in the same plane the Gram determinant vanishes. This hypersurface will turn out to be useful later when determining boundary values for the differential equations.
The one loop integral family is defined as where a k ∈ Z and the propagators are given by See figure 10. The integral family has five master integrals. In the one loop case it is possible to find a basis, where all master integrals have uniformal transcendental weight.

JHEP04(2018)047
Such a basis is called UT or canonical basis. The master integrals g of an UT basis fulfill a particular nice differential equation [6] d g(x, s 1 , s 2 ) = ǫ dÃ(x, s 1 , s 2 ) g(x, s 1 , s 2 ) , (D. 5) where the ǫ dependence is completely factorized and the matrixÃ is a linear combination of logarithms with coefficients given by rational matrices. The set of all different logarithm appearing inÃ is called alphabet of the differential equation. A possible choice of a UT basis is given by 7 The alphabet consists of seven letters {log(1 + x), log(x), log(1 − x), log(s 1 ), log(s 2 ), log(s 1 x + s 2 ), log(s 2 x + s 1 )}. We remark that the last two letters appear as factors in the Gram determinant (D.2). We solve the differential equation in the variable x for arbitrary s 1 and s 2 in an ǫexpansion g = k=1 ǫ k g (k) using iterated integrals. Then the integration constants are functions of s 1 and s 2 . It turns out that they are determined by the analytic properties of the integrals in certain limits of x, and in terms of the two (trivial) bubble integrals In this way we get the full solution of the differential equation. The feature that boundary values can be obtained trivially from the differential equations and physical considerations appears to be rather general, and has been observed in many calculations, see e.g. [86]. In order to understand how to obtain the boundary values, we only need to consider two limits, as will be described presently. First we consider the limit of a straight Wilson line (v 1 = v 2 ), where we have x = 1. Physically we do not expect any divergent behavior near x = 1, but the letter log(1 − x) can in principle give rise to logarithmic divergences. Such divergences can be extracted from the iterated integrals using the shuffle algebra they fulfill. The condition that there are no such divergences then yields a linear system of equations. The second limit is approached when all external vectors lie in the same plane. Then the Gram determinant (D.2) vanishes, hence one factor (s 1 x + s 2 ) or (s 2 x + s 1 ) must vanish. We expect the integrals to remain finite in this limit. Assuming s 2 > s 1 > 0, it is sufficient to consider the limit x → −s 1 /s 2 . The calculation is then identical to the first limit, except that we leave the Euclidean region to analytically continue the iterated integrals. For this we extract the logarithms log(x) from the iterated integrals using the shuffle algebra and use log(−x + i0 + ) = log(x) + iπ. We need the one-loop master integrals to order O(ǫ 3 ) for the two loop renormalization. To this order all master integrals can be expressed in terms of harmonic polylogarithms with the arguments x, s 1 /s 2 , xs 1 /s 2 and xs 2 /s 1 .

D.2 Two-loop master integrals
We can express all planar Feynman diagrams needed in our calculation as subdiagrams of a single integral family. This is particularly straightforward to see when using dual coordinates. We define the integral family as G a 1 ,...,a 9 = e 2ǫ γ E d D k 1 iπ D/2 where a k ∈ Z and the propagators are given by (D.11) The subtopologies we need for the calculation of the master integrals appearing in the correlators are defined by G a 1 ,a 2 ,b 3 ,a 4 ,b 5 ,a 6 ,a 7 ,a 8 ,a 9 and G a 1 ,a 2 ,a 3 ,b 4 ,a 5 ,b 6 ,a 7 ,a 8 ,a 9 , where a k is an arbitrary integer and b k is zero or a positive integer. They are related to each other through interchanging v 2 ↔ −v 1 . In figure 11 both subtopologies are shown. Note that in the Feynman diagrams more subtopologies appear, but after integral reduction all master integrals can be mapped onto these two. For these subtopologies we have in total 47 master integrals.
As in the one-loop case, we derived differential equations for all of them. Fixing the boundary values is analogous to the one loop case. In order to fix all boundary values we computed the following integrals An interesting new feature of this calculation is that there is a sector of the differential equations that leads to elliptic polylogarithms. The relevant integral sector is shown in figure 12. The elliptic nature of the integral can be seen by considering the projection of the differential equations onto that sector, as we describe below. The latter contains two master integrals.
Our expectation was that the elliptic sector should be irrelevant for the computation of the divergent part of the correlator (4.14). Given this expectation, we made a choice of master integrals that have good infrared and ultraviolet properties, namely g 38 = ǫ 4 G 1,1,0,1,0,1,1,0,1 and g 39 = ǫ 4 G 1,1,−1,1,0,1,1,0,1 . Indeed, with this choice it turns out that the elliptic sector decouples from the integrals needed for our correlator, to the order in ǫ that was required. The remaining 37 master integrals with fewer propagators than the elliptic sector are in UT form (D.5), with the same alphabet as the one-loop master integrals. As in the one-loop case, up to the order in ǫ that is needed for the renormalization of the correlators (4.13), all master integrals at two loops can be expressed in terms of harmonic polylogarithms.
In order to bring this differential equation into canonical form, one needs to solve the differential equation at ǫ = 0. It is sufficient to find the solution for one of the two solutions, as the other one can be obtained via the coupled equations. Here we focus on the integral g 38 . This integral is both UV and IR finite.
The desired solution at ǫ = 0 can be found from the maximal cut of the integral [88,89], which is the natural generalization of the leading singularities [6] to the elliptic case. Using