The Sudakov form factor at four loops in maximal super Yang-Mills theory

The four-loop Sudakov form factor in maximal super Yang-Mills theory is analysed in detail. It is shown explicitly how to construct a basis of integrals that have a uniformly transcendental expansion in the dimensional regularisation parameter, further elucidating the number-theoretic properties of Feynman integrals. The physical form factor is expressed in this basis for arbitrary colour factor. In the nonplanar sector the required integrals are integrated numerically using a mix of sector-decomposition and Mellin-Barnes representation methods. Both the cusp as well as the collinear anomalous dimension are computed. The results show explicitly the violation of quadratic Casimir scaling at the four-loop order. A thorough analysis concerning the reliability of reported numerical uncertainties is carried out.


Introduction
Gauge theory is the language of the standard model of particle physics. Even more than 50 years after its first modern formulation by Yang and Mills [1], it remains a hard task to compute observables even in a perturbative expansion in the coupling constants. Beyond perturbation theory much less is known in general, with the particular exception of those theories that admit a dual description in the strongly coupled sector, such as that provided by the AdS/CFT correspondence [2]. This correspondence is by far best understood for the maximally supersymmetric, N = 4, Yang-Mills (SYM) theory based on the SU (N c ) gauge group, in 't Hooft's planar limit [3]. In this limit, where N c → ∞, remarkable simplifications occur. A lighthouse result in this direction is the Beisert-Eden-Staudacher equation [4]: this equation describes a certain observable known as the planar lightlike cusp anomalous dimension (CAD) at all values of the coupling in N = 4 and ties into integrability ideas. Weak and strong coupling expansions of this anomalous dimension have been matched to independently obtained results, see e.g. [5][6][7][8][9][10][11][12]. However, beyond the planar limit much less is known in general despite some very recent progress in [13,14]. For the cusp anomalous dimension no nonplanar correction had been computed in any theory until recently the first numerical result at four loops in N = 4 was presented by us in [15]. Beyond the AdS/CFT correspondence and especially at weak coupling, the N = 4 super-Yang Mills theory is also a time-tested sandbox to explore computational ideas, such as those motivated by Witten's twistor string theory [16]. These have ignited a long-running program to explore the space of on-shell observables, using on-shell methods. This article is a part of this program, aimed at computing the so-called Sudakov form factor in N = 4 SYM theory. This form factor can be used to isolate several interesting universal functions that are contained within it. Prime among these is the lightlike cusp anomalous dimension mentioned above. The cusp anomalous dimension plays a central role in the analysis of infrared (IR) divergences, as first pointed out in [17]. By extrapolating structures found through three loops a general conjecture was formulated in [18] that the nonplanar part of the CAD vanishes in any perturbative gauge theory. This became known as quadratic Casimir scaling of the CAD, see e.g. [18][19][20][21][22][23][24]. It was noted that the quadratic Casimir scaling may be violated to higher orders of perturbative expansion due to the appearance of higher Casimir operators of the gauge group [25], see also [26]. At strong coupling, this scaling is known to break down in N = 4 SYM [27]. In addition, instanton effects break the scaling [28]. Finally, ref. [15] disproved the conjecture in perturbation theory, see also the two recent works [29,30] which apply directly to quantum chromodynamics and also report violation of Casimir scaling.
The Sudakov form factor we consider is an observable which involves two on-shell massless states and a gauge invariant operator in the stress tensor multiplet in N = 4 SYM, (1.1) In N = 4 SYM, form factors were first studied thirty years ago in [31] and revived in the past few years at weak coupling  and at strong coupling [62][63][64]. There have been interesting recent studies of loop form factors of non-Bogomolnyi-Prasad-Sommerfield (BPS) operators [65][66][67][68][69][70][71][72][73]. For reviews, see the theses [74,75]. The present paper is aimed at elucidating the evaluation of the integrals that appear in the four-loop Sudakov form factor, with the expectation that the presented techniques can be applied more widely.
A key idea in this article is to make transparent the transcendentality properties of the Feynman integrals that make up the Sudakov form factor. It is known quite generally that at fixed orders in the expansion in the dimensional regularisation parameter ǫ of Feynman integrals only rational linear combinations of certain constants appear. These constants are known as multiple zeta values (MZV). In principle, also more general constants such as Euler sums can appear, but in the known terms of the Sudakov form factor through to three loops in N = 4 SYM, MZVs are sufficient. MZVs have a property known as transcendental weight which takes integer values. The number of independent MZVs is small for low weight, and a basis for these constants is formed by (see e.g. [76]) with increasing weight denoted by the subscripts. At fixed order in ǫ in a generic integral only terms up to a maximal weight appear. This maximal weight increases stepwise with the order of expansion. A special class of integrals is formed by those where only the maximal weight terms appear at each order in the ǫ-expansion. Assigning to ǫ a transcendental weight −1, these integrals have a well-defined overall transcendental weight, and will be referred to as uniformly transcendental (UT) integrals. The concept of transcendental weight is important as it is observed in many examples that in N = 4 SYM (and superstring theory) only terms with maximal weight appear. Although the origin of this is somewhat ill-understood, it at the very least makes for a useful tool. Moreover, a general conjecture [8,77] relates the maximal transcendental terms appearing in QCD directly to N = 4 SYM for certain quantities. An example of this kind is given by the quark and gluon form factors in QCD [78][79][80][81][82] and the Sudakov form factor in N = 4 SYM, where the maximal transcendentality principle was verified through to three loops and for all terms up to transcendental weight eight [37]. Examples for two-loop remainders were also found in [38,73]. For the three-loop form factor in N = 4 SYM, an expression in terms of UT integrals was obtained in [37]. In that case the master integrals were known analytically, facilitating the analysis. In the four-loop case generically the basis of UT integrals was unknown. In this article, it will be shown how to identify UT candidates systematically, and how to write the four-loop Sudakov form factor as a rational linear combination of UT candidate integrals. The result in the nonplanar sector will then be integrated numerically, yielding a large list of new integral results. What is surprising is the empirical observation that obtaining numerical results for UT integrals turns out to be substantially simpler than for generic non-UT integrals in the class under study, even though the integration techniques themselves do not make use of the UT property. The result is combined into the nonplanar cusp and collinear anomalous dimensions at four loops. The result for the cusp anomalous dimension was first announced by us in [15], while the result on the collinear anomalous dimension is new. We comment extensively on the numerics below, making use of the UT property to inform the error analysis. This article is structured as follows: section 2 contains a review and setup of the problem. In section 3, uniformly transcendental integrals are discussed both at the general level as well as for the specific observable under study. Of special interest is a general technique for obtaining candidate-UT integrals. The full form factor is expanded in terms of the UT basis in section 4. In section 5 we discuss the numerical integration of the appearing integrals in the nonplanar sector, present our results and perform a thorough analysis of the reported numerical uncertainties. We conclude in section 6. The article is supplemented by several appendices. In appendix A we give explicit results of the UT integrals in the nonplanar sector, while appendix B contains the parametrisation of the integral topologies in terms of loop and external momenta.
2 Review and setup 2.1 Infra-red divergent structure of the form factor in N = 4 SYM The perturbative expansion of the Sudakov form factor is fixed by supersymmetry and dimensional analysis as where p 1 , p 2 are two on-shell momenta, and q = (p 1 + p 2 ) is off-shell. In dimensional regularisation with D = 4 − 2ǫ, F (l) is a purely numerical function of gauge group invariants and ǫ. The coupling constant is normalised as We consider explicitly the SU (N c ) gauge group, although our results apply to any Lie group: up to the order considered, there is a one-to-one map from N c to Casimir invariants, see below.
The form factor is free of ultraviolet (UV) divergences, since the operator O in the stress tensor multiplet is protected. On the other hand, there are IR divergences due to soft and collinear singularities from the massless states. Setting q 2 = −1 and defining the normalised form factor as F = 1 + ∞ l=1 g 2l F (l) , the IR structure is described in the following form [83] (exponentiation structure of Sudakov form factor in more general theories was original studied in [84][85][86][87] where the leading singularity is determined by the cusp anomalous dimension (CAD) γ cusp , and the sub-leading divergence is related to the so-called collinear anomalous dimension G coll . 1 Besides analysing the IR structure of the form factor, one also has to investigate its colour structure. For a classical Lie-group with Lie-algebra [T a , T b ] = if abc T c and structure constants f abc , the quadratic Casimir operators in the fundamental (F ) and adjoint (A) representation are defined via (see e.g. [88]) There are different conventions of defining cusp and collinear anomalous dimensions in the literature. In our convention, the cusp anomalous dimension γcusp = l γ (l) cusp g 2l is the same as the function f (g) in [4].
is the number of generators of SU (N c ). The colour structure of the form factor at l loops in N = 4 SYM theory where matter is always in the adjoint representation is simply (C A ) l up to l = 3. Starting from four loops, the quartic Casimir invariant arises in addition, and hence in SU (N c ) gauge theory one has, besides the planar (i.e. N l c leadingcolour) contribution a nonplanar (i.e. N l−2 c subleading-colour) correction. Starting from six loops, additional group invariants appear [42].
The planar form factor has leading divergence ∝ 1/ǫ 2l at l-loop order. To compute the CAD, this function needs to be expanded down to ǫ −2 at l loops, combined together with higher terms in the Laurent expansion in ǫ from lower-loop contributions. As mentioned above, the first nonplanar correction starts at four loops, due to the appearance of a quartic Casimir invariant. The nonplanar part of the four-loop form factor takes the following form In particular, it has only a double pole in ǫ since, upon taking the logarithm in (2.2), this piece cannot mix with any planar contribution from lower loops. We emphasise that individual integrals that contribute to F NP will typically have the full 1/ǫ 8 divergence. The cancellation of these higher-order poles in the final result therefore provides a very strong constraint on as well as a non-trivial consistency check of the computation.
The form factor exhibits a Laurent expansion in the dimensional regularisation parameter ǫ. In this expansion, each term is expected to be a rational-coefficient polynomial of Riemann Zeta values ζ n , or their multi-index generalizations, ζ n 1 ,n 2 ,... , known as multiple zeta values (MZVs) (see e.g. [76]). In principle, even more general objects such as Euler sums can appear. However, as mentioned earlier, any analytically known piece of the form factor does not go beyond MZVs. The MZVs have a transcendentality degree which is the sum of their indices, i n i . Also, the regularisation parameter ǫ is assigned transcendentality −1. In N = 4 SYM, the finite part of the form factor is expected to have (maximal) uniform transcendentality, which at l loops is 2l, and which suggests that the CAD at l loops is of uniform transcendental weight 2l − 2. Indeed, the planar CAD at four loops in N = 4 SYM has transcendentality six and was computed as [9,10,12] (log F )  Figure 1. Integral topologies that contribute only to the planar form factor at four loops.
We will provide strong evidence that also the nonplanar form factor and in particular the CAD are of uniform transcendentality at four loops. A numerical result of the planar fourloop collinear anomalous dimension G (4) coll,P was obtained in [89]. Recently, also the analytic value of this quantity was presented [90].

Integrand and integral relations
The full four-loop Sudakov form factor including the nonplanar part in N = 4 SYM was obtained as a linear combination of a number of four-loop integrals in [42] based on colourkinematics duality [91,92]. Similar five-loop result was also obtained recently in [57]. For more details on colour-kinematics duality, see e.g. the lecture [93]. The explicit form of the integrals for the problem at hand can be found in [42]. There are 34 distinct cubic integral topologies, each with 12 internal lines, that contribute to the four-loop form factor. They are labelled (1) - (34) in [42] and we provide them in figures 1 -3 for convenience and further reference throughout the present paper.
The four-loop integrals take the generic form as where D i are twelve propagators and N (l i , p j ) are dimension-four numerators in terms of Lorentz products of the four independent loop and two independent external on-shell momenta. For each topology, one needs to pick six additional propagators (i.e. six irreducible numerators) to form a complete basis, and we label them D k , k = 13, . . . , 18. Such a choice is   i , where the superscript (n) indicates the topology, and the subscript i, i = 1, ..., 18 refers to the basis given explicitly in appendix B (see also Appendix C of [48]). We define D (n) 19 = (p 1 + p 2 ) 2 . Any given numerator can then be represented uniquely in the chosen basis.
A fundamental property of Feynman integrals, as those in equation (2.7), is that they obey integration-by-parts (IBP) identities [94,95], which follow from Working out the left-hand side gives a linear relation between different integrals. By solving linear systems of such equations, a generic Feynman integral can be expressed in terms of a set of basis integrals. This procedure is known as IBP reduction, and the set of basis integrals is also known as the set of master integrals. The form factor was expressed in terms of a set of master integrals in [48] using the Reduze code [96]. 2 The master integrals, however, have evaded full integration so far due to their overwhelming complexity. In addition, the full IBP reduction generically leads to coefficients that contain higher-order poles in ǫ. This requires to evaluate the master integrals to higher orders in the ǫ expansion, which further increases the size of the problem. In this paper a different strategy will be used by expanding the form factor in terms of a set of integrals which are each simple enough to integrate and have ǫ-independent prefactors. A particular subset of the IBP relations turns out to be very useful for our purpose. These are the IBP relations in which the coefficients in front of integrals are pure rational numbers and independent of ǫ. These 'rational IBP' relations have been obtained in [105] for the form factor presently under study as a subset of the full reduction. An example is shown in Fig. 4. Note that integral relations derived from graph symmetries are a particular subset of the rational IBP relations.

Uniformly transcendental basis
A key idea of the present study is to expand the form factor in a set of integrals that all have uniform transcendentality (UT), which will be referred to as UT integrals. Such a representation of the form factor will make manifest the expected maximal transcendentality property of N = 4 SYM, and has been achieved at three loops in [37]. As will be shown in the next section, the UT integrals turn out to be much simpler to integrate numerically compared to generic non-UT integrals of similar complexity, which is crucial for the computation at hand.
We will now turn to the question how to find UT integrals prior to explicitly computing them. There are basically three ways to show whether an integral is UT.
• A UT integral can be written in the so-called dLog form [106,107].
• The leading singularities, or equivalently, the residues at all poles of a UT integral must always be a constant [107][108][109]. This is conjectured to be a necessary and sufficient condition.
• A set of UT integral basis can lead to certain simple differential equations [110].
The last point regarding differential equations is not directly applicable to the Sudakov form factor at hand since it is a single-scale problem, and thus not 'differentiable'. See however [109,111] for a work-around by deforming an on-shell leg to be massive, thus creating a two-scale problem. Below we illustrate the first two UT properties using a simple one-loop example. Then their application to four-loop form factor integrals will be discussed.

Warm up: a one-loop example
A one-loop UT example is given by the following scalar triangle integral: This is a UT integral as evidenced by the explicit result in the ǫ expansion An important interesting point is that, despite that the integral requires regularisation to be well defined, the UT property can be understood in exactly four dimensions at the integrand level. In the following, we consider only the integral in four dimensions as It is convenient to parametrise the loop momenta such that only scalar integration parameters remain. The four-dimensional loop momentum can be parametrised as where p i = λ iλi , i = 1, 2 are the external on-shell momenta, and q 1 , q 2 can be chosen as the two complex solutions to for example, q 1 = λ 1λ2 , q 2 = λ 2λ1 . The integral in the parametric form is .
This can be written in the following dLog form 7) which, in terms of momenta, is equivalent to the form (3.8) The existence of the dLog representation implies that the integral is UT. As mentioned above, an alternative way to prove UT property is to consider the leading singularity. In the parametric form like (3.6), this is equivalent to check the residues at all poles of the integral: the constant leading singularity property translates to the simple pole condition for all parameters. Let us explain this in more detail. To check the simple pole condition, one needs to pick up a certain order of the parameters to take the residue. Consider the one-loop example, we can first take residue for α 1 at the pole of the first propagator . (3.9) Next, we take the residue for α 2 at pole 0 Further residue at pole α 2 = 0 → dα 3 dα 4 1 α 3 α 4 . (3.10) The remaining parameters obviously have only simple poles and the final residue is a constant. One needs to check all different orders of taking residues, and in all occurring poles. In any intermediate step, after taking a residue in a particular parameter, if one encounters other than a simple pole in a remaining parameter, the integral is not UT.
We would like to emphasise that the simple pole requirement should also apply to poles at infinity. To be more concrete, consider following simple examples. For the integral there is a double pole for α 2 at 1, thus it is not UT. As for another example, both have a double pole for α 2 at infinity, so they are not UT either. 3 The condition that only simple poles are allowed is related to the required existence of a dLog form where only a logarithmic singularity is allowed. However, it does not require to find the explicit transformation to the dLog form which can be very complicated in general. This simple-pole condition provides an essential constraint for the construction of UT integrals below. A related strategy was also used in [107][108][109].

Systematic construction
The aforementioned condition of simple poles is used here both to construct and as well as to check UT integral candidates. Given a four-loop form factor integral in four dimensions, there are 16 integration parameters, so in principle there are 16! ∼ 2 × 10 13 different orders in which the residues can be taken. Practically therefore, the simple-pole condition is verified by choosing a large number of random orders of taking residues. A non-UT integral typically fails the UT test well within a few hundred of such random checks.
This UT test strategy can be used to constrain the space of potential UT integrals when combined with an Ansatz for the numerator. For the four-loop form factor integrals, one can start with a linear Ansatz of mass dimension four numerators of a given topology. We then perform the above described residue tests. The requirement of absence of higher order poles provides linear constraints on the set of coefficients in the Ansatz by computing the residues of the higher order poles. Solving these linear constraints then yields a smaller Ansatz, and the process is repeated. For speed, it is better to first identify a sequence of residues leading to a higher pole by choosing the Ansatz coefficients to be random integer numbers. This sequence can then be used to derive the analytic constraint on the full Ansatz. Below we provide more technical details.
A full four-loop topology contains 12 lines (i.e. propagators). For many topologies one can simply ask the following question: which sets of 10-and 11-line integrals can be added to a given topology in the four-loop form factor such that the sum is UT? Suppose such a linear combination exists. Then it is obvious this is likely not unique: adding any linear combination of 10-or 11-line UT integrals will satisfy the same constraint. To find a basis for all these UT integrals, we take the form factor numerators N ff as they appear in the N = 4 theory as input and add to this the set of all 10-and 11-line integrals (of which there are 162). In this case the initial Ansatz looks like where D 19 := q 2 and the D's are the propagators as given in appendix B. Inserting the parametrisation (3.4) for each of the four-loop momenta gives a rational expression of 16 αtype parameters. Now one needs to identify a sequence of residues yielding a double or higher pole in the α parameters. Demanding that the pole becomes a simple one yields at least one constraint equation for the 162 + 1 parameters {a 0 , a j,k , b j,k }. Explicitly solving these linear constraint equations gives a smaller Ansatz. Now one repeats by again trying to find a sequence of residues that will yield additional constraints. After a number of iterations for the integrals in the case at hand, one has obtained a set consisting of one integral containing the 12-line parts of the form factor contribution and other integrals which contain at most 11 lines. These are a set of UT candidate (UTC) integrals. One can also ask the question which UT integrals with unit exponent propagators exist in a given topology, worrying later about expressing the N = 4 form factor in terms of these.
To answer this question one chooses a more general initial Ansatz such as (3.14) Here the simple pole condition will provide a set of linear equations of 190 parameters {b j,k }. The end-result for this wider initial ansatz will be a set which contains all possible UT candidate integrals in a given topology (with unit exponents for the propagators). If, after deriving constraints with a certain number of random checks and no new further constraints are found in typically a few hundred more random pole checks, the remaining Ansatz contains a set of good UT candidates.
The choice of initial Ansatz is dictated to a large part by practical ease of subsequent numerical integration. For many public codes, the numerator of integrals is in general preferred to be a product of two factors, each quadratic in momenta. If a single such integral is to contain the full 12-line parts of a particular integral topology k, a necessary but not sufficient condition is to check that the irreducible numerators of a given integral form a product form separately. Concretely, one sets all propagators of this topology to zero, and verifies if a product form emerges for the irreducible numerators. In our chosen set of expressions, the propagators of a topology are always the first 12 entries (see Appendix B), so to check is: where the product form is a quadratic function of D i , i = 13, ..., 19. This condition is satisfied for all topologies in the four-loop form factor under study, except for topologies (12), (17), (19), and (26). Note this condition is independent of the exact choice of propagator basis. If this condition is satisfied, then the smaller Ansatz approach of form factor integral plus 10and 11-liners has a chance of sufficing. This is usually much quicker and more transparent.
If the 12-line parts do not have a product form, the larger Ansatz must be used. Examples of both possibilities are, for instance, topology (19) which does not have a product form, and topology (23) which does. From a generic set of UT candidates UTC i , the product form can be found by solving the following equation for non-trivial parameters λ, α and β which are rational numbers. This is a quadratic set of equations, obtained by matching coefficients of products of D's. Since we are interested in integrals that can be used to express the form factor in, more constraints can be added to the problem for specific purposes. For instance, the constraint can be added that the twelve-line parts match known form factor numerator contribution in the topology under study, Note this constraint only makes sense in a topology where the form factor has a product form on the left hand side of equation (3.15). Alternatively, one can simply demand one specific coefficient to be unity, This in particular avoids finding trivial solutions to the general problem in equation (3.18) (λ i = α j = β k = 0). This constraint is particularly useful when looking for very general solutions to the quadratic problem, matching only to some terms appearing in the form factor. Finally, one can add manifest graph symmetry constraints on the UT candidates: this we did in almost all cases. Which constraint to use in a particular situation depends on the generality of the solution sought for.
Having set up the quadratic problem (3.18), the first step is to solve the linear subproblem for λ. Then, one can impose graph symmetry patterns on the product form. The remaining set of quadratic equations can be analysed completely, or a particular solution can be guessed by computer algebra. 4 In several cases, it can be shown that no solution to a given problem exists. In these cases, after exhausting all options, one can widen the Ansatz in equation (3.18) by adding a linear combination of ten-line integrals (which are expected to be simple to integrate). These cases can be clearly seen in the results in section 4, e.g. (4.7) -(4.9). Also, sometimes residual parameter-containing solutions to the product-form problem are obtained. In these cases educated guesses were employed, aimed at as parametrically simple as possible integrals.
The result is a list of product-form UT candidates for each topology. The ones listed in this article have all individually been checked to pass at least 10, 000 simple residue checks, giving ample evidence for their uniform transcendentality. As will be discussed later, checking a set of found integrals individually also serves as a useful cross-check on computational errors.

dLog forms
Writing a four-loop integral in dLog form will give a direct proof of UT property. However, the construction of a dLog form for a generic four-loop form factor integral is a difficult task, and hence this method is more suitable to show the UT property of a given integral rather than to derive a UT numerator.
A useful strategy to construct a dLog form is loop by loop [107,108]. With proper numerators, all one-loop triangle and box integrals can be written explicitly in dLog forms. For example, the three-mass box is known to have a dLog form (see e.g. [107], k 1 is massless, K 2 and K 4 are massive) with given numerator which is the Jacobian of the quadruple cut of the box, such that the leading singularity is a kinematics-independent constant. So when there is a three-mass sub-box in the four-loop integral, one can write this sub-box in a dLog form, and the remaining integral is a three-loop q p 1 p 2 Figure 5. Topologies for which it is straightforward to construct a dLog form.
integral involving a new propagator 1/N 3m . In some topologies, such a procedure can be done recursively loop by loop, so that the full integral can be written explicitly in the dLog form. This normally happens when the topology involves at least one box with at least one massless leg, and has some ladder structure. 5 Such cases include topology (1), (6), (13), (21), (23), (28), as shown in figure 5, whose dLog numerators are given, respectively, by

Full form factor in UT basis
Finding an expansion of the full form factor in terms of generic UT candidate integrals can be obtained by relatively straightforward linear algebra techniques. In addition, we discussed above how to find product-form numerators for candidate UT integrals. Combining the two involves quite a wealth of choices that can be made in intermediate steps. For the nonplanar form factor, we first found a linear combination of 12-line UT candidates which satisfies the requirement that the difference to the full result contained at most 11-line integrals.
Combining the remaining expression into UT candidates in the nonplanar sector was then a relatively easy task. In the planar sector, it turned out that more work was required. An obscuring factor is the existence of many relations between different integrals from the rational IBP relations. A choice that works is given below. This choice was driven by the attempt to find as simple expressions as possible and to express the end-result in as small a number of integrals as possible. This includes both aiming at graph-symmetric expressions as well as trying to find an expansion involving only small integer or half-integer expansion coefficients. This necessarily involves some heuristics. It would be very interesting to find concise target integral expressions more easily, ideally driven by integration convenience or accuracy, but this would lead us beyond the scope of this work. One important result that follows is that both the planar as well as the nonplanar sector of the form factor can be expressed in terms of rational (i.e. ǫ-independent) linear combinations of UT integral candidates. We regard this as strong evidence for the maximal transcendentality of the form factor. By extension, this implies maximal transcendentality for the cusp and collinear anomalous dimensions at the four-loop order in maximal SYM theory, both in the planar and nonplanar sectors. Moreover, the smallness of the expansion coefficients clearly suggests this expansion is natural. In the nonplanar sector we have checked explicitly that the form factor integrals found originally in [42] when taken as complete topologies can only be expressed in terms of UT integrals in one unique combination of the 14 topologies: the one in which they appear. This provides a cross-check on the symmetry and colour factors.

UT integrals for the nonplanar form factor
Below we list 23 UT integrals I (n) 1 − 23 that combine into the nonplanar form factor. The superscript (n) denotes the twelve propagators from topology (n) in Fig. 2. In this notation, we only have to list the numerator of each integral. Moreover, each integral I (n i ) i gets multiplied by a rational pre-factor c i according to (4.1) The nonplanar form factor is then obtained as where the prefactor 48/N 2 c = 2×24/N 2 c is the normalisation stemming from the permutational sum of external legs and the colour factor [42], and the UT integrals are We note that integrals I 1 − 11 , I 12 − 18 , and I 19 − 23 , are 12-, 11-, and 10-line integrals, respectively. The integral I in topology (25) is the only one which does not carry the symmetry of the topology explicitly. This was done to arrive at a simpler form to integrate. In general topologies (25) and (26) are the hardest topologies to find UT integrals which are reasonably compact. Note that topologies (31) through (34) do not appear: there are no UT candidate integrals at all in these topologies.

UT integrals for the planar form factor
Similar to the nonplanar part, we also provide an expansion of the planar form factor in terms of 32 UT integrals I  The planar form factor is then obtained as where the prefactor 2 is the normalisation stemming from the permutational sum, 6 and the UT integrals are (as in nonplanar case, we only indicate the numerator) 19 )D (4.32) 3 + D 10 − D 18 − D 19 )D 3 − D 9 + D 19 )(−D 3 + D 5 + D 6 + D 17 + D 19 ) (4.34) 13 D 17 + D 18 + D 19 )(D − D 11 + D 14 + D 18 + D 19 ) − D 16 + D 19 )(−D 17 + D 19 ) 15 + D 17 ) 6 Note that unlike the nonplanar case, there is no color factor contribution.
In the present work, we choose a numerical approach. While numerical integration of the four-loop form factor integrals remains quite hard for generic numerators, we make the surprising empirical observation that UT integrals are numerically much easier to integrate than generic numerators of the class under study. We may offer an intuitive explanation for this. The constant leading transcendentality criterion used to find candidate UT integrals guarantees that these integrals have very mild singularity properties. An algorithm like sector decomposition is bound to be more efficient in cases where internal singularities are simpler. Note however that sector decomposition algorithm works in Feynman parameter space, whereas constant leading singularity criteria are applied in parametric form like (3.6). Whatever the precise origin, the relative simplicity of UT integrals is a boon for explicit computation, leading to a remarkable reduction in intermediate expression sizes and integration times. Moreover, the obtained coefficients in the expansion appear to be numerically much smaller than for generic integrals; this is beneficial for reducing potential cancellation errors.
Because of the physical motivation, we will only focus on integration of the integrals in the nonplanar sector of the form factor. We leave integration of the integrals in the colourplanar sector to future work, mostly because all terms of the latter through to O(ǫ −1 ) are already known: The ǫ −{8,6,5,4,3} poles are dictated by contributions from lower loops according to eq. (2.2), and the cusp [4,9,12] and collinear [89,90] anomalous dimensions are already known analytically.

Mellin-Barnes representations
Mellin-Barnes (MB) representations constitute a powerful tool for evaluating Feynman integrals [116][117][118]. They rely on the fact that one can factorise sums of terms at the cost of introducing line integrals in the complex plane. The basic formula reads The curves are usually straight lines parallel to the imaginary axis whose constant real parts are chosen such as to separate all left from all right poles of Γ-functions. This is achieved by choosing the real parts of all MB variables w i together with that of ǫ such that the arguments of all Γ-functions have positive real parts. The poles in ǫ are then extracted by analytical continuation to ǫ → 0, for which several algorithms exist [118][119][120][121][122]. Subsequently, the terms can be Laurent-expanded about ǫ = 0 and integrated, which proceeds mostly numerically with MB.m [121], but also examples of analytical evaluation of MB integrals exist [123,124].
To get from a loop integral to an MB representation, one first introduces Feynman parameters, e.g. like After integration over k 1 the remaining terms can be factorised using eq. (5.1), and subsequently integrated over the x i via x i ) = Γ(a 1 ) Γ(a 2 ) · · · Γ(a n ) Γ(a 1 + a 2 + . . . + a n ) .
The procedure is then repeated until all loop momenta are integrated out. In our case where no kinematic thresholds are present one has to obtain positive definite terms at all stages of the calculation if q 2 is space-like (we put q 2 = −1 for definiteness). Moreover, all terms in the ǫ expansion are real. For planar topologies this so-called loop-by-loop approach is always applicable, and we will refer to MB representations coming exclusively from positive definite terms as valid MB representations. However, valid MB representations for a given loop integral are not unique, even their dimensionality can differ depending on the order the loop momenta are integrated over.
For crossed topologies, the situation is more complicated as one encounters cases in which the loop-by-loop approach yields polynomials in the Feynman parameters x i which are not positive definite, even in the absence of kinematic thresholds. Consequently, the MB integrals will be highly oscillating and hence their numerical evaluation will be difficult to handle, although steps in this direction have been undertaken [125][126][127][128].
One way of circumventing this problem in the case of crossed topologies is to not integrate over the loop momenta one by one, but to simultaneously integrate over all loop momenta. This is done by means of the Symanzik graph polynomials U and F [129][130][131][132], which at L loops are homogeneous of order L and L + 1, respectively. In the absence of kinematic thresholds they are positive definite and hence automatically lead to valid MB representations. The price to pay is the fact that the number of terms in U and F scales as L!. As L increases this therefore quickly leads to MB representations that are too high-dimensional to be integrated in practice. A partial remedy to this problem is to group the lengthy sum of terms x i x j . . . in U and F into a short sum of products (x i + . . .)(x j + . . .) . . ., see our example below.
To take advantage of both the loop-by-loop and the FU approach we apply a hybrid of the two approaches here (see also [125]): To keep the dimension of the MB representation at a manageable level, we first integrate via the loop-by-loop method over as many (say, ℓ) loop momenta as possible such as to not generate non-positive definite expressions. Afterwards, we use the FU method for the remaining 4 − ℓ loops. This ensures that we deal with positive definite terms at all stages of the calculation and hence automatically obtain valid MB representations, and still keep the dimension of the resulting MB integral moderate since the number of terms in U and F only scales as (4 − ℓ)! instead of 4!. To give examples we look at different integral topologies from the nonplanar part of the four-loop form factor.
In topology (24) (see Fig. 2) for instance, we can integrate out the box that is attached to the external p 2 -line. Contrary to expectations the loop-by-loop approach fails if one tries to integrate out next the box attached to the external p 1 -line or any other loop. Consequently, the remaining three loop momenta have to be treated by means of the FU approach. The obtained MB representation is 21-fold and hence at the edge of what is doable in practice. We did some checks through to O(ǫ −4 ), where at most 7-fold integrals appear in the numerical evaluation with MB.m [121]. Topology (23) with numerator (ℓ 3 − p 1 ) 2 2 is even worse since we didn't find a single loop that can be integrated over before one is enforced to switch to the FU method. Consequently, the MB method was not applied to this topology. In topology (25), on the other hand, one can integrate over the two boxes attached to the external p 1 and p 2 -lines, respectively, and switch to the FU method afterwards. Still, the obtained MB representation is 20-fold. We use it for some checks through to O(ǫ −5 ), where at most 5-fold integrals appear in the numerical evaluation.
One example where the hybrid method works particularly well is topology (30) with numerator (ℓ 3 11 . Let us therefore give more details on this case. After shifting the loop momenta according to . We now integrate out the two boxes parameterised by the loop momenta k 2 and k 1 . Using eqs. (5.1) -(5.3) this introduces eight MB parameters and we are left with the propagators that are raised to various powers which can also depend on the MB variables. At this stage it becomes obvious that all terms in the numerator except k 2 4 can be treated as inverse propagators, giving rise to shifted propagator powers. While in principle also k 2 4 could be treated in this way, it would introduce an additional propagator and hence longer F and U polynomials. We therefore use the formulas in section 3.2.4 of [133] (which are based on ideas in [120,134]) for explicit numerator factors in case of k 2 4 . The F and U graph polynomials for this two-loop topology can now be written down, and a crucial step consists of writing the expanded sum of terms as a short sum of products which still happens to be positive definite. One obtains and hence U and F are sums of three and five terms only, respectively. Moreover, they do have various factors such as (x 4 + x 5 + x 6 ) and (x 3 + x 7 ) in common. The factorisation of the terms in F and U via eq. (5.1) now proceeds in several steps. The final integration over the x i is performed by means of eq. (5.3) and requires the introduction of an additional regulator δ in order to avoid a Γ(0) in the denominator. We choose to add δ to the power of x 4 . After application of Barnes' lemmas the resulting MB representation is 14-fold. The subsequent analytic continuation δ → 0 is done with MB.m [121] and the dimension is reduced to 13. This integrand is attached to the arXiv submission of the present article.
The package MB.m is also used for most of the remaining steps: Analytic continuation to ǫ = 0, expansion in ǫ, application of Barnes' lemmas, and numerical integration. After these steps, at most six-fold MB integrands appear through to O(ǫ −2 ). At O(ǫ −1 ) the maximum dimension of the integrand is seven. We use the algorithms CUHRE and VEGAS from the CUBA library with up to 2.1 billion sampling points. The result is given in appendix A.

Sector decomposition
Sector decomposition [132,135] regularises the ǫ-expansion of Feynman integrals by performing a blow-up at singularities of the Feynman parameter representation of a given integral. Sector decomposition has been implemented in several public codes, e.g. FIESTA [134,[136][137][138] and SecDec [133,139,140]. For our production runs, we have used the FIESTA code, with cross-checks in simpler cases from SecDec. After resolving the singularities, a list of remaining integrals is obtained. These could be integrated analytically in principle, but most often these are integrated numerically. The numerical integration is performed using mainly the VEGAS algorithm [141] as implemented in the CUBA library [142], with some crosschecks using CUHRE and DIVONNE from the same library. For our production runs, we In the course of computation several tricks were used to speed up computation and to control the arising errors. The sector decomposition programs involve choices of how to regularise the integrals, which are encapsulated in different strategies for resolving the singularities. This is a feature of which the problem at hand benefits a lot since the occurring integrals are complicated multivariate expressions. Whenever it finishes, FIESTA's "strategy X" typically leads to smallest sector counts which we will take to be a proxy for the ease of integration. In cases where this strategy fails for one or more sectors, one can split the computation into those sectors treated with strategy X and a remainder tackled with "strategy S". This can be done using the option "SectorCoefficients" in FIESTA. A further trick to use is that of graph symmetries. These can be used to gather exponents of several sectors into a single one, with a numerical pre-factor counting the number of sectors related to the base sector. For choosing the base sector, one first runs FIESTA on all sectors, selecting representatives which have the smallest sector count as a proxy for simplicity.
Note it is very important to verify that the integral in question has the explicit graph symmetry used; otherwise a wrong result may be the consequence. Here the UT properties of the integrals offer some protection: if an error with respect to graph symmetries is made during computation, the obtained final result is typically not UT and this manifests itself for instance by a non-vanishing ǫ −7 coefficient. Moreover, in these cases the numerical value of the coefficients tends to grow very fast with increasing orders of ǫ. In addition, if a graph symmetry is misused one cannot rationalise the coefficients of the ǫ expansion as described below.

Nonplanar cusp and collinear anomalous dimensions
We gather the numerical results for all integrals needed for the nonplanar part of the Sudakov form factor in appendix A. When combined to give the Sudakov form factor, the results are gathered in table 1. Errors are added in quadrature, see below for the rationale behind this. Due to the high precision of the computation at order ǫ −8 , there is no sensible reported error in FIESTA. Note that in table 1 the prefactor 48/N 2 c in (4.2) is not included.
As mentioned above, physics dictates that the coefficients of orders ǫ {−8,−7,−6,−5,−4,−3} vanish in the final result, which is numerically indeed the case and provides a strong consistency check of our computation. The coefficients of order ǫ −7 must even vanish in each of the 23 UT integrals separately. The orders ǫ {−8,−6,−5,−4,−3} are in most cases non-zero in individual integrals but cancel in the final result. As described below, the precision of the orders ǫ {−8,−6,−5,−4} is good enough to translate the reported numbers into small rational multiples of {1, ζ 2 , ζ 3 , ζ 4 }. After doing so, these orders also vanish analytically in the final result of the nonplanar form factor.
As can be seen from table 1, the first non-zero term is at order ǫ −2 . The result 1.60 ± 0.19 has a statistical significance to deviate from zero of 8.4σ. Adding individual uncertainties linearly to account for potential systematic effects would yield 1.60 ± 0.58; still significantly non-zero. 7 We will argue below that there is no evidence for systematically underestimated error bars in our calculation.
Translating the result of the order ǫ −2 of the nonplanar form factor into a result for the sought-after nonplanar four-loop CAD yields for gauge group SU (N c ) where the prefactor 3072 = 2 × 24 × 64 is the normalisation stemming from the permutational sum, the colour factor [42], and the denominator of (2.5), respectively. Compared to the planar result γ (4) cusp,P = −1752ζ 6 − 64ζ 2 3 ∼ −1875, we observe that the nonplanar CAD has the same sign. If we use N c = 3, its value becomes γ where the prefactor 384 = 2 × 24 × 8 has the similar origin as γ (4) cusp, NP above. Interestingly, compared to the four-loop planar collinear AD result, G coll, P ∼ −1240 [89,90], we observe that the nonplanar central value result +(6904 ± 1248)/N 2 c indicates the sign is different; it is also different from the sign of the nonplanar cusp AD above. This is a new feature comparing to all known planar results in which collinear AD always has same sign as cusp AD. 8 Note that our result is in tension with a vanishing result at the 5.5σ level. The largest contribution to the error budget within the integrals at this order comes from I (27) 8 , which contributes ∼ 1.86, followed by four integrals which contribute between 0.95 and 1 each, whereas all others are below 0.75. We mention that the linearly summed error is obtained as −17.98 ± 11.89.
To improve the reported uncertainties significantly within our numerical approach would come at a high price, both with respect to computing time and power, since the resources required for pushing this computation through sector decomposition are fairly large. It would certainly be interesting though to confirm the sign of the collinear AD.

Rationalisation
Since the used integrals pass all applied UT checks, their ǫ-expansion is expected to be UT. Assuming that MZVs are sufficient and no genuine Euler sums occur, for the orders ǫ {−8,−6,−5,−4} it is expected that the numerical coefficients can be written as a rational number times {1, ζ 2 , ζ 3 , ζ 4 }. Hence, by dividing the numerical result by the appropriate MZV constant, a numerical result is obtained which should be expressible as a rational number. For the case at hand, we typically have at least five to six digits available and the found integers have on the order of three digits in numerator and denominator. This indicates that the obtained rational numbers are reasonable, which gets supported by the fact that their contribution in the final result of the nonplanar form factor cancels exactly. In appendix A the results of the rationalisation are listed.
For ǫ {−3,−2} the UT property still holds, but at these orders there are two MZVs of transcendentality 5 and 6 respectively. For weight 5 these could for instance be taken to be ζ 2 ζ 3 and ζ 5 , and one can attempt a solution with the PSLQ algorithm [143], for instance through Mathematica's command FindIntegerNullVector. The appropriate integer relation then contains three unknowns: one for the numerical result, and two for the MZVs. For integer coefficients to be reliably isolated one needs much more digits in these cases, certainly more than 10. Since we have typically only four to five digits available at these orders, the PSLQ algorithm is currently not feasible. Moreover, many of our numerical results were obtained using sector decomposition where the price of integration roughly scales quadratically with increasing precision. This makes PSLQ unfeasible for the coefficients at orders ǫ {−3,−2,−1} within the numerical setup employed here. It would be highly interesting to obtain high precision numerics at these orders, or even better of course analytic results that do not rely on PSLQ.

Error analysis
Since numerical integration methods are used, a thorough discussion of the errors in these integrals is called for. Both for MB as well as for sector decomposition methods an error is reported. As is well-known, if an efficient MB representation can be found, the error in its integration is in general small, especially compared to sector decomposition. For the integrals at hand typically a difference in precision of three to five digits arises. Hence, the discussion here will focus on sector decomposition.
FIESTA employs the CUBA [142] integration library. Although we have cross-checked some simple integrals as well as leading expansion coefficients of more complicated ones, most of the coefficients needed for the cusp anomalous dimension at order ǫ −2 were obtained using exclusively the VEGAS [141] algorithm. VEGAS employs an adaptive sampling algorithm.
It should be noted that the integrals under study do not have any physical singularities, and do not have to be analytically continued, two common sources of error. For sufficiently many evaluation points, the VEGAS error is of Gaussian type. To check that this regime is reached, one evaluates the integrals for several evaluation points settings. In the Gaussian regime, the error scales as 1/ √ eval points. For all integrals in the set integrated here, this was reached very quickly. In rare cases involving much more complicated integrals, it has been reported in [144] that the error in FIESTA can be underestimated. In those cases the central value of certain coefficients changed outside the reported error with increasing evaluation points. We have checked for this as well, and have never observed variations outside of reported error upon increasing the number of evaluation points for the integrals under study. Several simpler integrals have been computed using SecDec with the DIVONNE and CUHRE algorithms as a further crosscheck. More cross-checks for integral I and I 16 follow from available MB results, as well as an exact result for integral I Finally, physics provides a strong cross-check of the numerics. The leading coefficient of the nonplanar form factor should be of order ǫ −2 , while individual integrals generically contribute from order ǫ −8 . Hence, in the sum there should be numerical cancellations between the integrals to give zero within error bars for the first six orders of expansion, down to ǫ −3 . With the errors added in quadrature and the result for the sum of the central value, one can compare to the exact answer, 0, for these coefficients. These results are contained in table 1 and clearly indicate that reported errors are not underestimated, giving further support for our error analysis.
In total, the above analysis shows that the errors reported by FIESTA are stable and in general conservatively estimate the errors for the form factor integrals in the present study. This strongly indicates that the final error for CAD is not underestimated either, and hence there is no need to manually inflate the reported uncertainty. Conservatively, we will interpret the FIESTA reported error as the standard deviation of a Gaussian error. For a true single We can see that all ratios are larger than unity, which suggests that the FIESTA errors are conservative estimates. Besides, we find that the deviation of FIESTA results from PSLQ results are both positive and negative, which indicates that there is no source of systematic errors. standard deviation in a Gaussian error, one would expect deviations from the true result to exceed the standard deviation of the Gaussian distribution roughly 32% of the time, while here this never occurs. As a consequence of the error interpretation, the obtained errors are added in quadrature. For reference, also the result of adding errors linearly is provided, which is recommended in cases which involve a small systematic error. However, we emphasise that there is no sign of systematic errors in the case at hand.

Discussion and conclusion
In this article a set of tools and techniques have been discussed for the integration of four-loop form factor integrals, especially focussed on the nonplanar sector of the Sudakov form factor in maximally supersymmetric Yang-Mills theory. This sector contains among others information on the nonplanar correction to the cusp anomalous dimension. Four loops is the first time a nonplanar correction enters into the form factor as well as into the cusp and collinear anomalous dimensions. Although conjectures existed that the CAD vanished generically in gauge theories, our results, first announced in [15], show this is not the case. In this article we also present the first numerical result for the nonplanar collinear anomalous dimension. The numerics of especially the latter result leave quite some room for improvement. Even more interesting would be to obtain an analytic result. Besides settling conjectures, of much wider interest is how the results reported in this article were obtained: the tools and techniques are certainly applicable to a wider context than just this particular computation in this particular theory.
Inspired by similar computations in the literature [106][107][108][109], an algorithm was presented to find complete sets of uniformly transcendental integrals in a given set of topologies. The algorithm is based on the conjecture that these integrals always have constant leading singularities. Importantly, the algorithm stabilises to a result in finite time in our current Mathematica implementation. A surprising amount of uniformly transcendental integrals were found for each integral topology for the problem at hand. With some combination techniques, a set of integrals was obtained to express the maximally supersymmetric form factor in. However, the number of UT integrals involved in this physical problem is much smaller than the total number of UT integrals in each topology. This points towards applications of these integrals beyond maximal supersymmetric Yang-Mills. Intriguingly, the numbers obtained are comparable to the total number of IBP master integrals. It would be very interesting to explore this further, but this will have to involve IBP-reducing the pure, non-supersymmetric Yang-Mills form factor, which is beyond currently (publicly) available technology.
Having obtained a suitable basis of UT master integrals to express the form factor in, the next step is the integration of these integrals. A pleasant surprise is the observation that even though many integration techniques such as sector decomposition spoil UT properties in intermediate steps, the UT integrals appear to be much easier to integrate than generic integrals in the form factor class. Within sector decomposition, this manifests itself in term counts which are an order of magnitude better. This in turn leads to much more compact expressions in the integration steps which lead to much improved performance in both speed and accuracy. Intuitively, this corresponds well to the notion that UT integrals are inherently simple. More mathematically, the absence of higher order singularities in the integrand in parametric form (as discussed in section 3.1) translates very likely to less singular integrands in Feynman parameter form. This in turn should then explain the observed much improved behaviour of sector decomposition methods. It would be interesting to explore this further, especially a criterion which would allow one to decide if an integral is UT in Feynman parameter form would be highly desired. Since there are considerably fewer integrations in Feynman parameter form than in parametric form, this is potentially even much more powerful.
Special attention is paid to the numerical integration of the form factor integrals in the nonplanar sector. Apart from the central value, the error analysis in numerical applications is important. Here the UT property of the integrals informs the error analysis. The integration of leading coefficients allows one to check the error analysis by using the PSLQ algorithm to find the exact value of the integrals. This combination of number theory and numerical integration shows that the errors reported by FIESTA are in general very conservative estimates. Added to knowledge of a single exact integrals and several results obtained using Mellin-Barnes integrals, this gives comprehensive evidence for our error analysis for the computation of the nonplanar cusp and collinear anomalous dimensions at four loops.
the Collaborative Research Center 676 "Particles, Strings and the Early Universe". GY is supported in part by the Chinese Academy of Sciences (CAS) Hundred-Talent Program, by the Key Research Program of Frontier Sciences of CAS, and by Project 11647601 supported by National Natural Science Foundation of China.

A.1 UT integrals with 12 lines
For the UT integrals we use the parametrizaton in terms of loop momenta from [48] and the normalisation used by FIESTA, i.e. we work in D = 4 − 2ǫ-dimensional Minkowskian spacetime and our integration measure is e ǫγ E d D ℓ/(iπ D/2 ) per loop. Moreover, we set (p 1 + p 2 ) 2 = −1 and suppress the fact that the ǫ-expansion continues in all equations. Below we give our numerical results as well as the PSLQ up to ǫ −4 order. is known analytically from [109]. Our numerical results obtained by MB and FIESTA agree with the analytical one well within error bars.

Topology 21
Topology 23 This result was obtained with MB. FIESTA performs poorly in this topology. (A.23)

B Basis of propagators and numerators
This appendix contains the basis of 12 propagators and 6 irreducible numerators, which are used in section 4.2. The numbering of the equations corresponds to the topologies in figure 1 -2. In each case, the first twelve entries parametrise the twelve propagators of the respective integral and the last six entries the chosen numerators. We have defined q = p 1 + p 2 .