The full angle-dependence of the four-loop cusp anomalous dimension in QED

The angle-dependent cusp anomalous dimension governs divergences coming from soft gluon exchanges between heavy particles, such as top quarks. We focus on the matter-dependent contributions and compute the first truly non-planar terms. They appear at four loops and are proportional to a quartic Casimir operator in color space. Specializing our general gauge theory result to U(1), we obtain the full QED four-loop angle-dependent cusp anomalous dimension. While more complicated functions appear at intermediate steps, the analytic answer depends only on multiple polylogarithms with singularities at fourth roots of unity. It can be written in terms of four rational structures, and contains functions of up to maximal transcendental weight seven. Despite this complexity, we find that numerically the answer is tantalizingly close to the appropriately rescaled one-loop formula, over most of the kinematic range. We take several limits of our analytic result, which serves as a check and allows us to obtain new, power-suppressed terms. In the anti-parallel lines limit, which corresponds to production of two massive particles at threshold, we find that the subleading power correction vanishes. Finally, we compute the quartic Casimir contribution for scalars in the loop. Taking into account a supersymmetric decomposition, we derive the first non-planar corrections to the quark anti-quark potential in maximally supersymmetric gauge theory.


INTRODUCTION
Low-energy gluons lead to large contributions in gauge theory scattering amplitudes and cross sections. The ability to predict and resum these enhanced terms is critical for phenomenology. It has been understood that scattering amplitudes factorize into a process-dependent finite part and a universal divergent part [1]. The latter is determined by a set of anomalous dimensions and by the soft anomalous dimension matrix, which describes the renormalization of products of Wilson lines pointing in the particles' directions [2]. A deeper understanding of infrared divergences is critical for developing manifestly finite approaches to cross sections [3], or for deriving quantitative predictions for the finite part of scattering amplitudes from symmetries of the S-matrix [4].
Many studies have focused on the planar limit of massless gauge theories, where significant simplifications occur [5,6], and where one may glean insights from integrability [7]. However, non-planar terms in QCD may be numerically dominant, and finite mass effects are important in many situations, for example when producing a top quark pair near threshold.
The multi-parton soft anomalous dimension matrix is known to two loops [8] and exhibits a simple dependence on the scattering angles. Whether or not this continues at higher loops is an open question. We will provide further insights into this issue by exploring the two-line case, which is given by the cusp anomalous dimension.
The cusp anomalous dimension depends on the Euclidean angle φ formed by the two Wilson lines. It is known to three loops [9][10][11]. Thanks to non-Abelian exponentiation, it depends on a reduced set of color factors, and truly non-planar corrections appear for the first time at four loops. They multiply quartic Casimir operators, and have recently been computed for the case of massless particles [12][13][14][15][16], thereby settling earlier questions about Casimir scaling. In this Letter, we consider the massive case and determine the matter-dependent quartic Casimir term in QCD. This determines the complete four-loop cusp anomalous dimension in QED.
Our work sheds new light on an observed iterative pattern, which allows to predict matter-dependent terms from lower loop orders [10,11]. Previous work showed that at four loops this pattern does not hold for all color structures [17]. The term we computed is precisely the one that captures the deviations from the conjecture in QED. We compare it quantitatively to the conjectured one, beyond the previously known limits.
We also derive new results in several physically important limits. Interestingly, we find that the first powersuppressed terms in the quark-antiquark threshold limit vanish. Moreover, we determine the first non-planar correction to the three-loop quark anti-quark potential in N = 4 super Yang-Mills (sYM).

METHODOLOGY
We compute the vacuum expectation value of the offshell Wilson lines in momentum space, using the conventions of heavy quark effective theory (HQET) [11]. We consider the gauge group SU (N c ), see [18] for our conventions, and perform calculations in covariant gauge, which allows us to check gauge invariance of the results. We obtain results for QED by adjusting the QCD color factors according to Figure 1: Feynman diagrams contributing to the four-loop quartic Casimir term B in eq. (1). [19], and by replacing the strong coupling strength α s = g 2 Y M /(4π) by the fine structure constant α.
We are interested in the matter-dependent quartic Casimir term of the four-loop cusp anomalous dimension, where B(x) and C(x) are the functions we wish to determine, x = e iφ , and where the number of (light) fermions and scalars (canonically coupled to the gauge field) is denoted by n f and n s , respectively. The dots stand for other color structures [18] which we do not consider. The Feynman diagrams relevant for the calculation of B(x) are shown in Fig. 1. We evaluate them, and the ones needed for C(x), in dimensional regularization, with D = 4 − 2 . The cusp anomalous dimension is obtained from the coefficient of the 1/ pole [18]. We relate the integrals appearing in the Feynman diagrams to a conveniently chosen integral basis [20], and then determine the basis integrals by the method of differential equations (DEs). In the process, we need to handle large systems of linear equations representing integration by parts identities, which we generate and solve using the codes FIRE6 [21] and LiteRed [22]. For most of the calculations we use FIRE6's finite field methods, and reconstruct the full x and D dependence using the techniques of [23]. We reconstruct only the DEs and the Feynman diagrams we need, as opposed to creating integral tables. Thanks to choosing an improved basis, relatively few finite field evaluations are required for this step.
The DEs for the integral families shown in Fig. 1 involve of the order of 500 integrals (needed to describe all integrals sharing the same or fewer propagators), with coupled sub-systems of size up to 17, and denominators of up to degree-20 polynomials in x and D. We wish to find an integral basis that significantly simplifies the differential equations [20]. To solve this complicated problem in an automated way we develop further the algorithm of [24] that relies on only a partial knowledge of a canonical basis. To provide input for the algorithm, we use several ideas: for up to nine propagators, we find candidate uniform weight HQET integrals by a position-space analysis [11], while for integrals with ten or more propagators we perform a leading singularity analysis with the help of the algorithm [14]. Moreover, we adapt [24] to scan over a larger set of candidate integrals automatically, and to use knowledge of all basis integrals that are known to be uniform weight. Finally, we also make use of the complementary method [25,26]. Given the information found by this means, we first transform the diagonal blocks (this is equivalent to considering integrals where propagators are replaced by delta functions) of the DE into canonical form. In a second step, we transform the off-diagonal blocks of the full DEs.
In this way we obtain a canonical form of the DEs for family (a), (b), (d) and (f), where f is the vector of basis integrals for each integral family, m k are matrices with constant entries, and (2) in terms of multiple polylogarithms [27] to the order in needed. The boundary values are taken from [17], or are determined from physical consistency conditions.
For the families (c) and (e), we proceed in the same way for all integrals up to nine propagators, and for all diagonal blocks, except for the ten-propagator sector shown in Fig. 2. For this sector we suspect that algebraic basis transformations are insufficient to bring them into canonical form (see Appendix). A canonical form can likely be achieved by enlarging the function space [27], but for the scope of this paper we can proceed in a simpler way.
We find that the issue can be bypassed in the computation of the cusp anomalous dimension, which only sees the 1/ pole of each HQET diagram. This suggests that complicated terms in individual basis integrals might drop out of the final answer. Having this in mind, we organize the basis in a way such that each basis element only needs to be computed to O( 0 ). Remarkably, the new basis forms a closed differential system up to O( ) corrections (which are irrelevant for Γ cusp ), where the matrix a is finite as → 0, and it is nilpotent. To the order needed, the solution is given by multiple polylogarithms.

Full result for the QED cusp anomalous dimension
We obtain the following result for the four-loop cusp anomalous dimension in QED, where and is the one-loop function. The first term in eq. (3) comes from propagator-type diagrams [28], while the term B, our main new result, comes from the Feynman diagrams shown in Fig. 1. We write it in the following way that makes the rational dependence on x manifest. The functions B i are given by linear combinations (in Q) of multiple polylogarithms of transcendental weight four to seven. We provide them as computer-readable ancillary files. Remarkably, only four of the integration kernels of eq. (2) appear, namely α = {x, 1 ± x, 1 + x 2 }.
Previously it was conjectured [10] that B(x) equals This was found to be inconsistent with the small angle expansion, but curiously it approximately agrees numerically with the exact answer [17]. We can now evaluate our exact result B(x) − B c (x) for any value of x, using [29]. Remarkably, we find that the shape of B(x) is very well described by B c (x). It turns out that one can improve the quantitative agreement by adjusting the proportionality constant in eq. (7). Fig. 3(a) and Fig. 3(c) show that −0.458A(x) approximates B(x) within 6% in a large part of the kinematic regions for Minkowskian angles ϕ = −iφ, and Euclidean angles, respectively, and similarly for the imaginary part in the region x ∈ (−1, 0), as shown in Fig. 3(b). These last two regions can be interpreted as below and above threshold for the production of two massive quarks, see e.g. [30]. For the real part for x ∈ (−1, 0), the deviation can reach 25% near x = −0.2, see Fig. 3(d). Still, we find it remarkable that the simple one-loop function A(x) captures the main features of the four-loop result.
Let us now expand our novel results in several interesting limits. This serves both as a check and allows us to produce novel results, by computing additional terms.
2. In the light-like limit x → 0, we find The first log x term agrees with the light-like cusp anomalous dimension [12,13], and the finite part is new.
3. In the anti-parallel lines limit x → −1, the cusp anomalous dimension is related to the quark-antiquark This relation is true up to beta function terms [11], which do not affect the quartic Casimir term considered here. We find where l 2 = log(2). The pole term is in agreement with [34]. It is remarkable that the subleading terms in eq. (9) start at O(δ) only.

Quark-antiquark potential in N = 4 sYM
The quark-antiquark potential in N = 4 sYM is known up to two loops for both the bosonic and supersymmetric Wilson loop [35]. At three loops, it has the following color dependence, We can determine V 2 for bosonic static charges by using a supersymmetric decomposition, as in [14]. Taking into account the known results for the gluon and the fermion quartic Casimir terms [34], we only need the scalar contribution. We obtain the latter from our result C(x) in eq. (1). (The full formula for C(x) is provided in an ancillary file.) We find This constitutes the first non-planar correction to the quark-antiquark potential in N = 4 sYM. It would be interesting if integrability methods [36] could be extended to this case.

DISCUSSION AND CONCLUSIONS
We computed the matter-dependent quartic Casimir term of the four-loop angle-dependent cusp anomalous dimension, in a generic gauge theory involving massless matter fields. In particular, this determines the full QED result for this quantity, and includes for the first time contributions from light-by-light scattering diagrams.
Our calculation revealed new structures compared to previous three-loop results: the answer contains four as opposed to previously two rational structures. The function space is given by iterated integrals, and compared to three loops, there is one new integration kernel, namely d log(1 + x 2 ), in addition to d log x, d log(1 − x) and d log(1 + x). This relative simplicity is remarkable given that intermediate steps (and higher order terms in the dimensional regulator) contain further integration kernels, cf. eq. (2). This hints at better approaches that avoid the complicated intermediate terms.
The information on the function space is valuable input for bootstrap approaches [37]. In particular, our result implies constraints (via collinear limits) on the multi-leg soft anomalous dimension matrix for massive particles, which currently is known at two loops only [8].
We analyzed the novel four-loop result numerically, and found that, surprisingly, it is described within a few per cent by a rescaled one-loop function (and up to 25% for the real part in the above threshold region). Clearly this could be improved even further by using some of the known limiting behaviour, as in [38], for example. It would be interesting to understand why the approximation works so well here, in view of other problems where fully analytic results are not yet known, such as scattering processes with many mass scales.
To obtain all color contributions of the complete QCD cusp anomalous dimension at four loops, only two further contributions are needed [18]. The first one is the planar limit, which is conceptually easier compared to our calculation. The HQET integrals we computed should cover a large part of the integrals needed, and are available upon request from the authors. The second one is the gluonic quartic Casimir term. Thanks to our calculation of the scalar terms, the latter can equivalently be obtained from the non-planar N = 4 sYM contribution. This opens up novel approaches to this problem [14].

Integrals beyond multiple polylogarithms?
Two of the integral families contain the integral sector shown in Fig. 2. It contains 11 coupled basis integrals. Following the procedure outlined in the main text, we find the following obstacle to reaching a canonical form for the DE. At order 0 we encounter the following two-by-two system (after eliminating spurious singularities [25]), where y ≡ − (1−x) 2 x and y ± = 3 √ 3 e ± iπ 6 . While the residue matrices at the singular points 0, y ± have eigenvalues 0, at y = ∞ we find eigenvalues -1 and 1. We were not able to find a balancing transformation [25] that sets all eigenvalues to zero. It would be interesting to prove that this two-by-two system cannot be balanced via algebraic transformations, as it would mean that this integral sector involves functions beyond multiple polylogarithms.