Master integrals for the four-loop Sudakov form factor

The light-like cusp anomalous dimension is a universal function in the analysis of infrared divergences. In maximally ($\mathcal{N}=4$) supersymmetric Yang-Mills theory (SYM) in the planar limit, it is known, in principle, to all loop orders. The non-planar corrections are not known in any theory, with the first appearing at the four-loop order. The simplest quantity which contains this correction is the four-loop two-point form factor of the stress tensor multiplet. This form factor was largely obtained in integrand form in a previous work for $\mathcal{N}=4$ SYM, up to a free parameter. In this work, a reduction of the appearing integrals obtained by solving integration-by-parts (IBP) identities using a modified version of Reduze is reported. The form factor is shown to be independent of the remaining parameter at integrand level due to an intricate pattern of cancellations after IBP reduction. Moreover, two of the integral topologies vanish after reduction. The appearing master integrals are cross-checked using independent algebraic-geometry techniques explored in the Mint package. The latter results provide the basis of master integrals applicable to generic form factors, including those in Quantum Chromodynamics. Discrepancies between explicitly solving the IBP relations and the MINT approach are highlighted. Remaining bottlenecks to completing the computation of the four-loop non-planar cusp anomalous dimension in $\mathcal{N}=4$ SYM and beyond are identified.


Introduction
Maximally (N = 4) supersymmetric Yang-Mills theory (SYM) offers perhaps the best chance of truly solving an interacting four-dimensional quantum field theory. In addition, it is a proven stepping stone to pioneer new computational technology, as the large amount of supersymmetry renders many computations much simpler than their non-supersymmetric counterparts. In many cases, the resulting technology has transformed computational power in generic quantum field theories. For instance, in QCD, recent years have seen a boom of next-to-leading-order (NLO) computations considered unfeasible only a few years before, inspired ultimately by Witten's twistor string [1].
The focus of this paper is the computation of the Sudakov form factor, which is an observable that involves one gauge-invariant operator from the stress tensor multiplet and two on-shell massless states, In N = 4 SYM, this form factor was first discussed and computed at the two-loop order in [2]. It is noteworthy that the three-loop form factor was first computed in QCD [3], almost three years after the master integrals were found in [4]. This result was then adapted to provide the three-loop answer in N = 4 SYM in [5]. The latter computation cleanly shows that the computation within N = 4 SYM is technically much more straightforward, especially when employing modern unitarity-based methods instead of Feynman graph techniques.
Since there is only a single scale in the problem of computing the Sudakov form factor, the dependence on this scale is fixed by dimensional analysis. Hence, the form factor evaluates to a function which only depends on the dimensional regularization parameter ǫ (D = 4 − 2ǫ), the coupling constant and the number of colors N c . Moreover, the dependence on ǫ is governed to a large extent by the known structure of infrared (IR) divergences. The divergent structure through 1 ǫ 2 , for instance, is determined by the so-called cusp (or soft) anomalous dimension, which gets its name from its appearance in the computation of the light-like cusped Wilson line [6,7]. The cusp anomalous dimension is a function which is universal for a given quantum field theory.
In N = 4 SYM, an integral equation which determines the leading planar part of the cusp anomalous dimension in principle exactly was derived from integrability [8]. This is to date the most powerful and precise manifestation of the cluster of ideas known as the AdS/CFT correspondence [9]. However, much less is known about the non-planar part of the cusp anomalous dimension. At weak coupling, an immediate problem is that the first non-planar correction appears at four loops and has never been calculated to date, in any theory. Further motivation comes from a conjecture [10] that this correction may in fact vanish. This is based on the assumed completeness of the dipole contributions, which is consistent with the analysis of IR divergences in Yang-Mills theories explored through three-loop order, see also [11]. However, there is evidence from the Regge limit analysis in [12] that this naive dipole summation is not sufficient at four loops. It is important to settle this issue by explicit computation.
The integrand of the form factor was obtained in [13] using color-kinematic duality [14,15,16] as an ansatz generator. The coefficients in the ansatz were fixed by a choice of unitarity cuts, up to one remaining parameter. To evaluate the integral, a major next step is to perform the integral reduction, in particular using integration-by-parts (IBP) reduction [17,18]. This turns out to be highly non-trivial due to the complexity of four-loop integrals. The success of this reduction is a major achievement of this paper. In particular, this reduces the challenging non-planar cusp anomalous dimension problem down to the computation of an explicit set of master integrals.
Based on the IBP reduction result, surprisingly, the free parameter in the integrand turns out to be truly free: it drops out of the full form factor after the reduction. This completes the determination of the integrand of the N = 4 SYM form factor at four loops. In other words, we prove that there is a one-parameter family of four-loop form factor integrands which satisfy color-kinematic duality.
The master integrals are also explored independently by algebraic techniques. While the results are consistent with the IBP reduction in majority, some rare examples of mismatch present interesting open questions to understand further. The results from algebraic methods, with some qualifications, provide the master integrals of four-loop form factors in generic quantum field theory, including QCD.
The remaining problem of computing the four-loop form factor in N = 4 SYM is the evaluation of the master integrals. Many techniques are available to compute Feynman integrals, see [19] and references cited therein for an overview. However, the basis integrals of the four-loop form factor are rather non-trivial. They include non-planar four-loop integrals with up to quadratic irreducible numerators, the majority of which, to our knowledge, have never been integrated. We present partial numerical results and briefly survey the challenges involved in completing the computation, while leaving the full evaluation to future work.
This article is structured as follows. First, the connection between form factor and cusp anomalous dimension is reviewed in section 2. This section also introduces some of the language to study Feynman integrals and discusses the issue of graph symmetry. In section 3, we report on the IBP reduction of the Feynman integrals in the four-loop form factor in N = 4 SYM based on the Reduze package [20]. The independence of the sum on the free parameter reported in [13] is elucidated. We then proceed, in section 4, to study the basis of master integrals using algebraic techniques explored in [21]. A brief survey of the numerical evaluation is given in section 5. Appendices contain information on a three-loop cross-check and an alternative choice of basis integrals. In the supplementary material section, there are several ancillary files which contain results for the basis of master integrals.

Review and setup
In the following, extensive use will be made of the results in [13] for the four-loop form factor in N = 4 SYM, in particular of the graphics and tables in section 5 of that paper. The numbering of the integral topologies refers to that paper.
In the course of writing this paper, a typo was discovered in Table 5 of [13]: the entry for the color factor of integral topology 26 should be 24N 2 c δ a 1 a 2 . In particular, this integral topology does not contribute to the planar cusp anomalous dimension.

Relation between form factor and cusp anomalous dimension
The Sudakov form factor is an observable which here involves one gaugeinvariant operator from the stress tensor multiplet and two on-shell massless states, In the present paper, we focus on the four-loop Sudakov form factor in N = 4 SYM, for which the integrand was obtained in [13]. 1 1 In N = 4 SYM, form factors were first studied thirty years ago in [2] and revived in the past few years at weak coupling [5,22,23,24,25,26,27,28,29,30,31,32] and at strong coupling [33,34,35]. See also recent developments concerning form factors of non-Bogomolnyi-Prasad-Sommerfield (BPS) operators [36,37,38,39].
The Sudakov form factor for the stress tensor multiplet can be written as the tree-level form factor times a scalar function, It is this function F 2 which will be computed. In momentum space, F 2 is a function of the only scale in the problem, q 2 = (p 1 + p 2 ) 2 , the coupling constants, group theory Casimirs and ǫ. Here, q is the momentum associated with the gauge-invariant operator. Note that it does not obey a mass-shell condition, in contrast to the momenta of the gluons, p 2 1 = p 2 2 = 0. The cusp anomalous dimension γ cusp is related to the two-point form factor through the universal exponentiation of IR divergences [40,41,42,43], which, in dimensional regularization, takes the form The cusp anomalous dimension can be calculated as a perturbative series in the coupling constant. Displaying the dependence on group theory quantities up to fourth order, it reads where C A is the Casimir which, for the SU(N c ) gauge group, simply has value N c . The quantity d 44 is a particular group theory Casimir invariant which, for SU(N c ), has the value N 4 c + 36N 2 c . In previous works, the constants a 1 , a 2 , a 3 and a P 4 were calculated in perturbation theory [7,44,45,46,47] 2 and using integrability [8]. In 't Hooft's planar limit [51], N c → ∞ with g 2 N c fixed, it is clear that the first subleading color correction to the cusp anomalous dimension occurs at four loops.
From the general exponentiation formula in equation (4), this non-planar correction can be isolated from the non-planar part of the form-factor at four-loop order at order 1 ǫ 2 in the ǫ expansion. Since the integrals present in the form factor naively diverge as 1 ǫ 8 , this gives an important check on any computation of the non-planar form factor: the first six orders in the 1 ǫ expansion must vanish. The planar correction to the cusp anomalous dimension at four loops can be expressed as a sum of rational coefficients times either ζ(3) 2 or π 6 [45,46,47], which are both constants of transcendentality six. Based on the maximal-transcendentality principle [48], it is natural to conjecture that the four-loop non-planar correction a NP 4 is given by a combination of the same transcendental numbers, if not vanishing.

Feynman integrals
An L-loop Feynman integral with n so-called indices a 1 , . . . , a n is an integral which can be written as where D i are inverse propagators. For future reference, we define the absolute value of the sum over positive and negative indices to be t and s respectively. The parameter s is the numerator power and t is the propagator power. In this article, the indices are strictly integer valued. Integrals which have the same non-zero indices for propagators are said to belong to the same sector. In the special case where all propagator powers are one, the integral is referred to as the corner integral of the sector.
The problem under study in this article consists of a series of four-loop twelve-propagator integrals. For each integral topology, one can construct a basis by starting with a parametrization of the loop momenta of the graphs. Then, one needs to pick six additional propagators to get a complete set. These additional propagators are important for expressing irreducible numerators. Given a numerator, one can express this uniquely into the basis. The choice of numerators is by nature somewhat arbitrary. If one focuses solely on additional propagators which are products of differences of two momenta, then this can be done by trial and error. A more systematic approach is needed, if the aim is to include graph symmetry.
As an example, let us consider integral topology 26. Its propagators can be parametrized as: This topology has two independent permutation symmetries, as will be discussed later. It can be checked that the following six additional propagators parametrize all irreducible numerators: In section Appendix C, choices of propagators are explicitly given for all 34 topologies. A useful step toward computing the Feynman integrals is the so-called α representation, see e.g. [19]. The result is an integral over α parameters, one for each propagator, where the normalization constant, which is unimportant for our purposes, has been suppressed. U and F are certain polynomials of the α parameters of homogeneity L and L + 1, respectively. See [19] for further discussion and [52] for a Mathematica code to compute them from given propagators. As is customary, the labels on the α i parameters correspond to the position in a given ordered list of propagators, such as the one in equation (7) for integral topology 26. Equation (9) can cover negative indices, that is, numerators as well.

Graph symmetry
Several of the graphs in the set under study have a graph symmetry: the graphs are mapped to themselves under a permutation of some of the edges and external legs. Exploiting graph symmetry can be of great help in simplifying and cross-checking computations if it is present.
In [13], graph symmetry was indicated briefly by the symmetry factor. These symmetries can be checked or found by using the built-in capabilities of Mathematica. These permutations can include permutations of the external legs. An alternative strategy is to compute the U and F polynomials of the α parameter integral, equation (6). Then, one can simply check all possible permutations and see if they leave these polynomials invariant. Due to the factorial growth of this problem, the feasibility of this analysis strongly depends on the algorithm: a naive implementation ran out of steam for ten propagators. A faster algorithm is presented in [53].
A symmetry of the graph of an integral can correspond to a symmetry of a corresponding Feynman integral. For a given parametrization of the loop momenta of the graph, a permutation symmetry of the graph may be undone by a shift of the loop momenta, combined with a permutation of the external legs. If the latter leaves the integral invariant, the result is a symmetry of the Feynman integral.
For the example of topology 26, the two independent permutations are generated from as can be checked explicitly in equation (7). Since the integrals, by dimensional analysis, only depend on (p 1 + p 2 ) 2 , they are left invariant by the permutation of p 1 and p 2 .
The momentum map given above corresponds to a simple permutation pattern for a given list of propagators. Let us consider, for instance, the scalar integral in topology 26, with the propagators as listed in equation (7). The two independent permutation symmetries of the graph correspond to the following cycles: Together with the trivial permutation and their product, these form a fourelement representation of the permutation group. A given propagator can be mapped to other propagators in the set, but only if they are in the orbit of the given permutation. The orbits of the above permutation group are This can be used to simplify the sector decomposition method employed in FIESTA, for instance, as will be shown in section 5.1. Given a permutation symmetry of a particular topology, ideally, additional propagators would be found which yield a complete basis of propagators and also incorporate the permutation symmetry. One can search for symmetric numerators systematically. Permutation symmetries either leave propagators invariant or interchange two propagators. Hence, the six propagators to be added to the set contain either zero, one, two or three pairs of interchanges and, consequently, six, four, two or zero invariant propagators. Starting from a generic linear combination of all four loop momenta and two external momenta, one can construct the most general pair combination as well as the most general invariant propagator in terms of a number of free variables. Since propagators are only determined up to a sign by the appearance of the square, there are two times two possible (pairs of) polynomials which come out of this. Hence, all possible appearances of the permutation symmetry can be parametrized.
Then, one checks if the resulting set of eighteen propagators is linearly independent. In several cases (e.g. topology 20), no set exists at all. In addition, there are integrals with multiple permutation symmetries which do not permit a choice of six additional propagators with all permutation symmetries manifest. In cases in which there are sets of numerators with explicit permutation symmetries, typically many free variables exist. These may then be fixed by hand aiming for propagators as simple as possible. Care should be exercised not to choose parameters such that accidental linear dependencies in the total set of propagators are introduced.
For the example of topology 26, the momentum maps in equation (10) do not map the additional numerators in equation (8) into a permutation of themselves. The following set of additional numerators has explicit permutation symmetry: This example shows a general principle: the price of manifest permutation symmetry can be a much more complicated numerator structure. It would be very interesting to fix guiding principles for choosing sets of numerators.
The considerations above result in two sets of propagators. For each integral topology, there is a choice of eighteen propagators which can be used to express any integral in this particular topology. In particular, any numerator can be expressed in terms of the basis. In the first set, the aim was to include numerators as simple as possible, while, in the second, permutation symmetries were taken into account as much as possible. Obtaining the symmetric set is the result of a fairly lengthy computation. The simplest choice in the first set is listed in section (Appendix C). All results reported here were obtained with this set, unless explicitly mentioned otherwise.

IBP reduction and its output
The set of Feynman integrals in equation (6) is over-complete, as there are relations between different Feynman integrals. A particular example are the IBP relations [17,18], which follow from Working out the right-hand side leads to a linear relation of different Feynman integrals. By solving a system of such equations, one may express a general Feynman integral in terms of some basic integrals. A different way of phrasing the problem is to envision the system of IBP relations as a large matrix equation, with the integrals combined into a vector. The standard way of solving a problem of this type is Gaussian elimination. However, the output of this algorithm for a non-invertible matrix problem, such as the one under study here, depends on the ordering of the integrals in the vector. This is the essence of Laporta's algorithm for IBP reduction [54]: one picks an ordering of the integrals. This ordering should be such that complicated integrals are expressed in terms of simpler ones in general. For instance, integrals with smaller values of the parameter s (measuring numerator powers) are preferred, as are integrals with less numbers of different propagators. In addition, smaller values of the parameter t (measuring propagator powers) are to be preferred. Within these general choices, in practice, the exact criteria can differ between two different implementations of Laporta's algorithm. The results of these two implementations are, therefore, in general not the same, but related by a change of basis.
Given a system of equations, the output of Laporta's algorithm is a reduction of all integrals in a given set down to integrals which cannot be further reduced from the given equations using the given ordering. The leftover integrals are known as master integrals. These integrals depend on the given set of equations, although typically, if the given set of equations is large enough, the set of master integrals tends to converge.
Various private and public implementations of Laporta's algorithm exist, such as AIR [55], FIRE [56,57,58] and Reduze [20,59]. See LiteRed [60,61] for an alternative approach to IBP reduction. We explored FIRE, Reduze and LiteRed in some detail for the four-loop form factor problem. Only Reduze succeeded in solving the problem, probably due to its better parallelisation implementation.

Relevant implementation details
Reduze is designed to run in massive parallel mode on a cluster. As an input, it takes a family of integrals. This family is treated as an ordered set: lower integrals in the set are considered simpler. Integrals are reduced as far as possible to lower integrals. This cuts down on the computational workload considerably, especially in a large-scale problem such as the one under study here. As parameters, Reduze takes a range of propagator powers (the t variable) as well as a range of numerator powers (the s variable) and constructs out of these all possible IBP relations for a given sector. These are then solved sector by sector, starting from the simplest sector, which is the bottom sector of the first member of the integral family. The number of relations tends to grow very quickly along either the t or the s axis.
A problem we encountered in the public version of Reduze in parallel mode was that it tended to crash when handling many large files on the cluster file system during the identity generation stage. Due to this problem it is unfeasible to obtain an IBP reduction with the public version of Reduze, and a fix is needed. The problem was that a file was reported as not present in the file system, while a manual check thereafter did uncover this particular file. This is due to the internal file handling structure of Reduze, 3 where all processes are allowed to read from and write to disc. As a work-around, the program can be forced 4 to sleep for a time-out of fifty seconds, followed by a 500-second time-out if the file is still not found. This two-stage timeout resolved the crashing problem due to disk space usage. The loss of productivity for a few minutes is a small price to pay to prevent the total collapse of the computation.
During the running of the four-loop reduction, Reduze still crashed occasionally when handling integrals in sectors with many, typically twelve or eleven propagators. These integrals require a large memory usage each, which can easily overwhelm one of the nodes of the cluster. Sometimes this can be due to bad scheduling of processes, with all root processes which consume most memory running on a single node; this situation can somewhat be avoided at start-up by passing scheduling instructions to the messagepassing-interface (MPI) program. 5 By contrast, for lower propagator count, the limiting factor for speed is mostly CPU time. These two situations require opposite numbers of allowed parallel processes, which in turn requires occasional input by the user during running. It should be noted that, although Reduze works in parallel, there is a limit to the number of cores which are assigned to a single process while still increasing computational speed. This saturation is sector dependent, but certainly under a hundred cores.
The chosen IBP relations are in the range t ≤ 12 and 0 ≤ s ≤ 2, with some extensions to t = 13 for specific integral topologies. Since Reduze features a choice of used relations, there can be problems with the reduction at the outer edges of these relations: there can be unresolved integrals. A simple check of this is to observe the file size of the reductions in a sector and compare it to the average for other integral topologies in a sector with the same numbers of propagators. A sudden increase in file size tends to indicate unreduced integrals: the unreduced integrals appear with much more complicated prefactors than the others, typically involving polynomial ratios with polynomial orders an order of magnitude above the norm. This is more than just a nuisance. The abnormal prefactors tend to be highly divergent in the limit ǫ → 0, typically ∝ 1 ǫ 15 or worse. Since the form factor itself is expected to diverge at worst as 1 ǫ 8 , this implies massive cancellations among the integrals. This particular problem occurred in the basis of section Appendix C in integral topologies 20 and 26, for integrals with at most ten propagators. This problem was resolved using consistency conditions, see two paragraphs down.
We reduced using the set of numerators in section Appendix C. A secondary computation using numerators with explicit permutation symmetries failed; at a certain eleven-propagator integral, the MPI protocol ran into an issue with maximal size of a single message, which is set by the maximum size of an integer on a 32-bit machine, see [62] for an explanation of the problem and a possible resolution. This is a result of a choice of protocol and not related to the physical architecture used. In principle, this problem could be circumvented by rewriting the Reduze code to handle messages of this size.
Since the set of numerators with the simplest structures was used, integral relations from permutation symmetries are not explicit. This allowed a simple check on the results: take a graph symmetry relation and use the IBP reduction. If Reduze found the relation, the result must be zero. We checked permutation symmetry on maximal-propagator integrals with upto-quadratic numerators. Surprisingly, one single additional relation was found which involved exactly all master integrals with abnormal prefactors. Plugging this relation into the reduction eliminated the abnormal prefactors, reducing them to much more normal-looking ones.

Results
The computation led to a result after several months of computation on a fairly large cluster. 6 The complete list of the master integrals obtained with this reduction is attached, both in the planar and non-planar sectors, as well as a sample result for topology 25. Overall, the integrals which appear after reduction are simpler. There are two integrals with two numerator powers. One of these contains only eight propagators and is easy to integrate numerically. The other, however, has the full twelve propagators and occurs in integral topology 26. It turns out that this particular integral can be integrated numerically (see section 5.1).
Any basis of master integrals listed is arbitrary to some extent. For a given sector, an obvious first choice for a master integral is its corner integral. Beyond this, there is a choice of a (typically single) numerator or doubledup propagator. The statistics in Table 1 relates to the direct output of Reduze, taking into account the one extra relation found from permutation symmetry. The statistics is split into those integrals which contribute to the planar and non-planar parts of the form factor.  In these tables, there is quite some double counting of integrals which appear both in the planar and non-planar sectors. Using the known result for the sum of the integrals in the planar sector, one can eliminate several of the integrals in the non-planar sector. This is more efficient than might be expected, as there are four integral topologies (21, 25, 30 and 31) which each appear with the same overall color factor, 2N 4 c + 24N 2 c . This can be used to eliminate these integral topologies from the non-planar sector in favor of the integrals in the purely planar sector and the known planar-form-factor result. The resulting sum contains more integrals (260 in total) than the non-planar sector by itself.
The choice of master integrals in the provided results is taken from the output of Reduze. It is possible to choose other sets of master integrals, in particular the ones with numerators, to aim at simplifying, for instance, the prefactors of the integrals in the form factor or to make use of a permutation symmetry of the scalar skeleton of the integral. In particular, one aim could be to make these as weakly divergent as possible in the limit ǫ → 0. The power of divergence in this limit determines the order to which the integral is to be expanded to get the cusp anomalous dimension.

Cross-check: Multiple reductions
The IBP reduction was run twice, with a second reduction involving a reshuffling of the integral topologies to put planar topologies first. It was checked that the two reductions commute: reducing with the second reduction followed by the first gives the same result as reducing with the first, after taking into account the extra relation found from permutation symmetries. Before taking into account this relation, there are seeming mismatches between several of the integrals which appear in the full fourloop form factor: these mismatches are all proportional to the extra relation. This is a consistency check for the IBP reduction.

Cancellations in the form factor after IBP reduction
The IBP reduction displays two separate patterns of interesting cancellations for the form factor result: one pattern takes place within the separate integral topologies, while the second is between different integral topologies.
The first pattern is that, for integral topologies 18 and 20, the contributions to the form factor vanish after IBP reduction. Generically, the Feynman integrals in these topologies do not vanish; there are twelve-propagator master integrals. However, the specific combination of the integrals in the N = 4 SYM form factor reduce to zero. Considering the complicated numerators for topologies 18 and 20, this is rather nontrivial and not obvious at all at the integrand level. It should be noted that these integrals already do not contribute to the form factor: their color factors are also zero. Such cancellations, however, do not happen for the other four zero-color-factor topologies, 8, 11, 15 and 16, whose integrals reduce to nontrivial results.
It would be very interesting to find an explanation for these vanishing integrals, as, in physics, a zero almost always has a physical and important symmetry explanation.
Another intriguing pattern emerges when tracking the dependence on the single free parameter in the form factor integrand which was left in [13]. There, the free parameter in the form factor multiplies a combination of 24 out of 34 different integral topologies which obeys all the physical constraints imposed in [13] separately. In particular, this combination is color-kinematic dual by itself. One may expect that, by applying other, more complicated cuts which are not exploited in [13], this parameter may be fixed uniquely. The surprise is that, when summed with the correct color and symmetry factors (and only then), the dependence on the free parameter drops out of the form factor after IBP reduction. This means that the colorkinematic-dual representation of the four-loop form factor contains a truly free parameter! This detailed cancellation involves an order of fifty different master integrals with up to ten propagators and occurs in both the planar and nonplanar sectors. Most of the cancellations of the coefficients of the master integrals involve integrals out of all of the eighteen integral topologies which involve the free parameter and have nontrivial color factors. As an intricate example, consider the coefficient of the master integral which occurs in topologies 12, 14, 15, 17 and 19. The coefficient in topology 12 reads: which is exactly −2 times the coefficient of this master integral in topology 19. Taking into account the symmetry factors of topologies 12 and 19, which are 2 and 1, respectively, this leads to a cancellation. Similarly, the coefficients of this master integral in topology 14 cancels against that in topology 17, while the color factor of topology 15 is zero. The other master integrals typically involve much more intricate cancellations.
It would be very interesting to reach a better understanding of the uncovered cancellations. On a practical level, this pattern completes the determination of the form factor integrand [13]. To our knowledge, this is the first example of a color-kinematic-dual integrand representation that contains a free parameter. It would be fascinating to see if there is any deeper physical meaning of the free parameter, for example, as due to some gauge symmetries. This free parameter can be set to any value at the start of the computation. A natural value is −1, which simplifies the naive appearance of the integrals. This cancellation provides strong evidence for consistency of both the integrand result in [13] and the IBP reduction performed in the present paper.

Counting master integrals using algebraic techniques
In [21], Lee and Pomeransky developed a method to count the number of master integrals using a technique adapted from algebraic geometry. This is partially based on previous work by Baikov [63]. This method only relies on the topology of the integral and, therefore, applies for arbitrary values of t and s. In particular, it is independent of obtaining an explicit IBP reduction and, hence, of the results of the previous section. This idea was also implemented in a Mathematica package called Mint [21]. In this section, we develop and apply this method to count the number of master integrals of four-loop form factors.
Similarly to the start of IBP reduction, one first needs to determine the unique physical sectors for a given integral topology. As described in [21], it is convenient to perform this step using LiteRed [60,61]. In particular, LiteRed determines the sector symmetries, which can be taken into account automatically in Mint. In the next step, the counting of master integrals is done independently for each physical sector.
Let us consider a physical sector of L-loop topology and m propagators. From the associated polynomials U and F in equation (9), one can define a new polynomial 7 which has polynomial degree L + 1. Using the duality of homology and cohomology (see [21] for more details), the number of master integrals is related to the number of independent integral contours. The latter depends on the analytic structure of the polynomial G and is characterized by the so-called proper critical points. These are the stationary points at which G does not vanish and are, thus, defined as the solutions of ∂G ∂α i = 0 (i = 1, . . . , m) and G = 0.
Proper critical points may or may not exist for a given G polynomial. If there is no critical point, this implies that the corner integral of the sector is reducible, so that this sector does not contribute a master integral. On the other hand, if critical points do exist, the corner integral is irreducible with respect to IBP. This condition of irreducibility was first shown in [63].
provided that the proper critical points are isolated. We will explain the meaning of isolated critical point shortly. In the case of non-isolated critical points, more work is required [21], see also an explicit example in section 4.1 below. Thus, the problem of counting master integrals is reduced to the one of computing the critical points of a system of nonlinear polynomial equations given by equation (18), which is usually highly nontrivial to solve. Fortunately, one can apply powerful algebraic approaches to make this job simple, based on the so-called Gröbner basis technique. Readers who are not familiar with the Gröbner basis and related concepts may consult [64] for a pedagogical introduction. Here, only the basic procedure is outlined.
The solution space of a system of polynomial equations is called affine variety associated with an ideal determined by the equations. For the system of equations (18), the ideal can be defined as where, in the last term, one introduces an additional parameter α 0 which forces the polynomial G to be non-vanishing at the critical point. The solutions of the affine variety can be obtained by computing the Gröbner basis of equation (20), gb(I) = {g 1 , g 2 , . . . , g k } .
Once the Gröbner basis is obtained, it becomes relatively trivial to find the solutions. In our problem, we only concern ourselves with the number of solutions. This, in practice, can be conveniently obtained by counting the number of irreducible monomials in the Gröbner basis. 8 The above-described procedure is implemented in Mint [21] and applies straightforwardly to many simple cases; for example, applying it to the three-loop form factor, we obtain the result summarized in Table 2. The numbers of master integrals are consistent with the reduction of LiteRed. Table 2: Master integral counting for the six three-loop form factor topologies given in Figure 5 of [13].

topology
(1) (2) (3) (4) (5) (6) # MIs via Mint 9 10 10 10 14 10 For the four-loop form factor, however, one encounters two problems when trying to apply the Mint package. The first problem is that the set of critical points can form an affine variety of dimension ≥ 1. This corresponds to the non-isolated critical-point case. Mint cannot deal with such cases automatically in its present version, but denotes them as Indeterminate.
The number of such cases for the integral topologies at hand are summarized in Table 3. Although the absolute number of these cases is low, they tend to occur in the more complicated sectors. They require some further work as described in [21]. We employed a similar procedure and present an explicit example in subsection 4.1. Another problem is that, in quite a lot of cases, Mint gets stuck at the step of Constructing Gröbner basis . . . . This is due to the complexity of four-loop integrals for which Mathematica cannot finish the computation of the corresponding Gröbner basis. There are many known other packages which can compute a Gröbner basis more efficiently. We applied Macaulay2 [65], which can be conveniently used in Mathematica through mathematicaM2.m [66], and we employed Singular [67] in several most complicated cases.
Based on Mint and with some extra effort to solve the above two problems, we obtained the numbers of master integrals for the 34 different topologies as summarized in Table 4.  (13) (14) (15) (16) (17) (18) (19) (20) (21) (22) (23)  Given the number of master integrals, Mint can also generate a set of master integrals. In most cases, a sector contains only one master integral, in which case the corner integral is chosen as the master integral. If there are more than one master integrals in a sector, Mint suggests integrals with higher powers of propagators (while Reduze typically chooses integrals with numerators). However, the choice of a double-propagator integral is rather heuristic; therefore, to find out whether it is a genuine master integral requires further checks. In our case, we could check them by using Reduze, as will be discussed shortly.
To obtain the master integrals for the four-loop form factor, we still need to combine the master integrals of Table 4 together. Obviously, same master integrals may appear in different topologies of Table 4. This may be detected by comparing their graph topology: if two graphs are isomorphic to each other, they are the same integral. After this step, we find that there are only 280 non-isomorphic master integrals: 244 with only simple propagators plus 36 containing a double propagator. 9 Here, one should be cautious that there could be relations between integrals which are apparently not due to graph isomorphism, see e.g. [20,53] for some explicit examples. Such relations can be found by considering Feynman parameter representations of integrals and checking if the two representations are related to each other by a permutation of Feynman parameters. A fast algorithm for such a check was proposed in [53]. Alternatively, one can also consider the so-called graph matroid technique to detect such hidden relations between different integrals, as implemented in Reduze [20]. For us, there is an easier way to detect such relations, namely applying Reduze to reduce the basis integrals obtained by Mint. We checked that all 244 corner integrals correspond to basis integrals as determined by Reduze. Thus, they are indeed master integrals.
As mentioned above, the choice of master integrals with double propagators by Mint requires further checks. We applied Reduze to the 32 doublepropagators integrals which do not appear in the twelve-propagator sector of topology 26 and found that four of the suggested master integrals are actually reducible. 10 This does not change the number of Mint master integrals, but simply means that a different integral from this sector must be chosen as a master integral. By changing the position of the double propagator, we indeed found irreducible integrals. The twelve-propagator sector of integral topology 26 is analyzed further below. Thus, we obtained the final list of 280 genuine master integrals.
We would like to emphasize that the counting of master integrals based on Mint only relies on the topologies of given integrals and applies to arbitrary numerators; therefore, it applies to any theory, including QCD. We remind that the counting of Reduze given in Table 1 concerns particular master integrals appearing in the reduction of the N = 4 SYM form factor result.

An interesting mismatch between Reduze and Mint
Some caution is needed when assessing the output of Mint: this package is known to under-report the number of master integrals in very rare cases.
To check this, we collected all master integrals appearing in the IBP reduction by Reduze, with t = 12 and s = 2 and with t = 13 and s = 1. It turns out that Reduze yields more master integrals. Many of them consist quarks, one likely needs all 280 master integrals. 10 Note that Reduze tends to choose integrals with numerators rather than double propagators. So, irreducible means that, in the reduction result of Reduze, there is an integral with the same propagators containing a numerator.
of twelve propagators and one to several numerators. They are at the edge of the reduction setting, and it is very likely that they may be reduced by including more IBP identities. These will not be discussed further here, details are available on request. What is unexpected is that there are nine corner integrals with only ten or less propagators, for example {0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 0, 0, 1, 0 with only eight propagators. Mint shows that they contain no proper critical points, while Reduze and FIRE5 take them as master integrals. Due to the low number of propagators, most related IBP identities should already be taken into account by Reduze. Therefore, it is likely that these integrals are truly master integrals. 11 This would mean that the counting method of Mint requires some modification, although it is correct in the vast majority of cases. We cannot rule out, however, that there are integral relations which render this set of nine corner integrals reducible. It would be very interesting to understand this issue and see which is right. In addition, there is one more master integral in the Reduze basis, which appears to be an artifact of the way we solved the permutation symmetry consistency condition.

Two worked-out examples
Here, we give further details by considering two examples which cannot be computed directly by Mint.
The quotient space of the last ideal is two-dimensional, while the others are all one-dimensional. So, naive counting would give seven master integrals. Topology 26, however, has permutation symmetry, see equation (11). It is not hard to see that the first four ideals are all related to each other by this symmetry; therefore, they contribute only one independent master integral. In total, one finds that there are four independent master integrals in this sector. For this integral topology, Reduze picks four twelve-propagator integrals as master integrals: the corner integral, two linear-numerator integrals and one with a quadratic numerator. Hence, the number of master integrals matches between Mint and Reduze.

Towards numerical integration
Given a set of master integrals, the next step is to perform the integration. A significant part of the master integrals obtained in the previous sections are very complicated and seem to be on or slightly over the very edge of what current integration methods can achieve. There are various methods to perform numerical integration in the literature, see [19] for an overview.
It should be said that, in principle, Mellin-Barnes (MB) representations, which can be used very efficiently in solving planar topology integrals [69,70,71,72], may be the optimal weapon of choice for four-loop integrals, if the difficulties in obtaining efficient MB representations for non-planar integrals without hidden singularities can be resolved. See [19], page 124, for a brief description of the problem. Progress in obtaining up to and including some three-loop non-planar representations in an automatic way was reported in [73]. This may be used, at least in principle, to obtain MB representations for those integral topologies containing either a bubble or a triangle with two on-shell legs.
We mostly explored numerical integration using sector decomposition [74]. In general, this is a method to turn a highly divergent integral into a sum over finite integrals, by resolving the singularities of the Feynman parameter integral. These finite integrals can then be integrated, for instance, using the CUBA library [75]. We are aware of three different implementations of this idea in public codes: sector decomposition [76], SecDec [77,78,79] and FIESTA [80,81,82]. Only the latter program in its latest versions (3.x) was able to handle some of the most difficult integrals of the problem at hand.

Relevant implementation details
There are various strategies for resolving the singularities of a given Feynman integral. FIESTA offers a choice of these. Most relevant here are strategies S [80] and X [74]. Strategy KU [83], which, in principle, is highly efficient in terms of sectors, did not resolve any of the sectors of the most complicated integrals (the corner integral of topology 25) in either of its three variants and was not pursued further. Strategy X, if it terminates, tends to yield better results in terms of the number of sectors than strategy S.
FIESTA splits the Feynman parameter integral in equation (9) into different primary sectors. In each primary sector, one particular Feynman parameter is smaller than the others. If the integrand is invariant under a permutation of the α parameters, then these primary sectors may be related. A particular permutation symmetry of the α parameters descends from the parent graph. This symmetry may then be used to reduce the computational workload. Let us recall the symmetry property of the scalar integral of topology 26 given in equation (11). From its orbits of the permutation group in equation (12), we can see that, instead of computing primary sectors 1, 2, 7 and 8 separately, one better just computes one of these, as the others will give the same numerical answer. In practice, it can matter which coefficient is picked, since the structure of the integrand for the α parameters is different; this can result in different numbers of sectors. The integration time varies; for integrals with twelve propagators, the preparation time alone is typically measured in days to weeks. By contrast, all integrals with less than eleven propagators are highly tractable using FIESTA, typically requiring less than a day.

Quadratic-numerator integral in topology 26
The quadratic-numerator integral in topology 26 is naively expected to be one of the hardest integrals in the full set. It can be integrated using FIESTA. Integral topology 26 has two permutation symmetries, see equation (10). To make use of this, one needs to pick a symmetric quadratic numerator. That is, one needs a numerator which is mapped to itself under both symmetries and which is related to the quadratic-numerator master integral in the output of Reduze by IBP reduction. Moreover, to facilitate the computation of the double derivative in the α parametrization, it is advantageous to pick a numerator that is linear in the loop momenta. It is easy to check that (l 3 · (p 1 − p 2 )) 2 (29) transforms into itself under the symmetries in equation (10). We checked that this is related to the quadratic master integral from the explicit IBP reduction. With the symmetries in place, one only needs to compute four out of twelve primary-sector coefficients in FIESTA. We picked the coefficients {4, 0, 4, 2, 2, 0, 0, 0, 0, 0, 0, 0} for the order of propagators in equation (7), basically choosing the lowest from the orbits in equation (28). The resulting integral can be integrated with some numerical effort, giving the result Preparing for the numerical integration takes about five days. The integration time itself strongly depends on the integration parameters in the Vegas integration algorithm of the CUBA library. 13 We checked by running with different numbers of evaluation points that the first five coefficients are numerically stable under variations of the relevant parameters, especially maxeval. In other words, it is likely that their errors are in the reported ranges.

Summary and discussion
In this paper, we obtained an explicit integral reduction for the fourloop Sudakov form factor in N = 4 SYM based on the integrand result in [13] and by exploiting the Reduze package [20]. This reduces the challenging four-loop problem, in particular the computation of the unknown non-planar cusp anomalous dimension, down to the evaluation of an explicit set of given master integrals.
We also obtained master integrals by using the algebraic method introduced in [21]. This provides a first non-trivial implementation of the method to a complete observable. Limits of the Mint package and their resolutions are explained. The resulting integrals are, with some qualifications, the master integrals of quite generic form factors. For instance, they are expected to cover all integrals which appear in the QCD result. By comparing with the results of IBP reduction, we found interesting mismatches, namely, existing masters integrals from IBP reduction are taken to be reducible by the Mint package. It would be important to understand the reason of this discrepancy.
Based on the reduction result, we furthermore observed surprising cancellations in the form factor. First, the form factor integrand obtained in [13] contains a free parameter. While it was unclear whether this parameter could be uniquely fixed by further physical constraints, we showed that the terms depending on this parameter all magically cancel upon using IBP identities. Therefore, this parameter is truly free and gives the first example of a color-kinematic-dual integrand with such a free parameter. Furthermore, two integrals, corresponding to topologies 18 and 20, turn out to be exactly zero, which is not obvious at all at the integrand level. It would be very interesting to have a physical understanding of these intriguing cancellations.
What remains to be done is to compute the master integrals. However, for many of them, this is currently just out of reach for publicly available (implementations of) integration methods. There is an important analogue to the three-loop case, which we would like to draw attention to. In that case, a reduction to master integrals was obtained in [4]. Their integration was first performed almost three years later in [3], using partly numerical methods. An analytic result was then reported in [84], with an important cross-check in [85]. Currently, the three-loop form factor integrals often appear as an example application of new integration and organization methods, e.g. in [86,87].
A new feature in the four-loop case is the appearance of several master integrals beyond corner integrals. These can be taken to be either doubledup propagators or single numerators. The latter can, in principle, be traded for higher-dimensional scalar integrals using the Feynman parametrization in equation (9). The problem is that these scalar integrals feature four dots: there is a sum over integrals with four doubled-up propagators. Without an IBP reduction of these integrals down to the basis, this is an unpractical result. An explicit IBP reduction with up to four dots would enable the use of several methods. For instance, one could use dimensional recurrences [88] to either compute the integrals numerically [89] or to express the form factor in terms of a quasi-finite basis [87]. An explicit quasi-finite integral reduction for the three-loop form factor has been announced [90]. The needed IBP reductions at four loops seem very hard to obtain within the current version of Reduze and simply out of reach for anything else on the public market. See [91] for a proposal of an improved algorithm which is particularly suitable for single-scale reduction problems.
There is a practical issue when calculating the form factor which concerns the choice of master integrals. There are several choices of master integrals one could make as for both linear numerators and doubled propagators. Which one to choose should then be guided by integrability, i.e. by which of these choices they can actually be integrated. Moreover, as for numerical evaluation, there are issues of precision and speed to be considered. One would expect doubled-up propagators to give simpler Feynman parameter integrals, while linear numerators tend to give less divergent coefficients. A canonical example of the two choices are the master integrals of the doublebox-integral topology. One can either choose a doubled-propagator basis element, as first done in [92]. Alternatively, one can use a linear-numerator integral, as in [93]. The latter has a less divergent prefactor as a function of ǫ. Which choice is the more practical must be decided on a case-by-case basis though: we observed cases where no difference in expansion order is present.
There are various other avenues to be explored which may lead to a full computation of the cusp anomalous dimension through the form factor route. One would be to explore MB representations of four-loop non-planar integrals. An obstacle here is that naive implementations of this idea quickly lead to high-order MB representations containing an order of one hundred MB integrations. Any algorithm which could curb this number would offer a way forward. While numerical evaluation may have the precision issue, analytic methods would be more desirable. Linearly reducible integrals may be evaluated analytically by HypInt [94]. This works particularly efficiently for integrals of less singular behavior [87]. An alternative route for computing the cusp anomalous dimension is to extend integration techniques based on differential equations [95] to form factors by making one extra leg off-shell and then taking the massless limit. Controlling this limit after integration might also be an obstacle. Apart from that, IBP reduction, a necessary step in constructing the differential equation, would be more challenging for deformed integrals. See [86] for an application to the three-loop case.
the form factor is which is to be compared to the known result, e.g. from [5], The difference between the exact answer and the numerical approximation is ∆F (A.3) Note that the final coefficient is only off by about one percent. With more numerical effort, this could be brought down more. For the purposes of this article, it is important that the calculational pipeline works, at least in principle.