{'fraction_non_alphanumeric': 0.03843730308758664, 'fraction_numerical': 0.014135685780298257, 'mean_word_length': 4.285999777950483, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 63, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multivariate Residues and Maximal Unitarity

We extend the maximal unitarity method to amplitude contributions whose cuts define multidimensional algebraic varieties. The technique is valid to all orders and is explicitly demonstrated at three loops in gauge theories with any number of fermions and scalars in the adjoint representation. Deca-cuts realized by replacement of real slice integration contours by higher-dimensional tori encircling the global poles are used to factorize the planar triple box onto a product of trees. We apply computational algebraic geometry and multivariate complex analysis to derive unique projectors for all master integral coefficients and obtain compact analytic formulae in terms of tree-level data.


Introduction
Potential discovery of new physics depends on our ability to compute precision cross section predictions for scattering of subatomic particles and in particular, a quantitative understanding of all relevant Standard Model processes which necessarily must be separated from the experimental data. Theoretical calculations carried out at the LHC start with treelevel amplitudes at leading order (LO) in perturbative Quantum Chromodynamics (QCD), whereas a combination of one-loop amplitudes and higher-multiplicity trees provide nextto-leading order corrections of quantitative reliability. Although computations of this type are very complicated, recent years have seen major advances at NLO, especially for processes with many final states. In the near future, theoretical calculations offered at NLO become insufficient to saturate the bound for experimental uncertainty. The upcoming frontier is therefore formed by NNLO computations and in particular, two-loop amplitudes, which are also relevant already at NLO for specific processes such as production of electroweak gauge bosons by gluon fusion, for example.
Historically, contributions to perturbative scattering amplitudes have been tracked pictorially by means of Feynman diagrams, which lead to precise mathematical expressions using the Feynman rules. Although this method gives an indispensable fundamental intuition for interactions between elementary particles, even simple problems beyond twoby-two gluon scattering beyond the tree-level become cumbersome due to the presence of a large redundancy in the theory needed in order to compensate for virtual intermediate states. In the last decade, strikingly powerful on-shell methods for amplitude computations at tree-and loop-level involving only physical information and analytic properties have emerged, the most prominent examples being the Britto-Cachazo-Feng-Witten (BCFW) [1,2] recursion relations and the generalized unitarity method, an enhanced version of the original unitarity method due to Bern, Dixon, and Kosower [4,5]. Their importance is reflected by the fact that all gauge theory and also gravity trees now may be constructed by just the Cauchy residue theorem and complex kinematics in three-point amplitudes, which then are recycled for loops.
Otherwise unattainable processes in massless QCD have been computed to at the twoloop order with four external particles using the unitarity method and other techniques including integration-by-parts identities [35][36][37][38][39][40][41]. In the last couple of years, new promising methods for two-loop amplitudes such as integrand-level reduction by multivariate polynomial division algorithms using Gröbner bases and classification of on-shell solutions by primary decomposition based on computational algebraic geometry have been demonstrated [53][54][55][56][57][58][59][60][61][62][63]. In particular, very recently the first results for five-gluon scattering at two loops in QCD were obtained by Badger, Frellesvig and one of the present authors along these lines [57]. The same technique has also been applied to maximal cuts of three-loop amplitudes [55].
In this paper, we will pursue amplitude computation at the level of integrated expressions in the framework of maximal unitarity, proposed by Kosower and Larsen [46]. Being an intense version of the generalized unitarity method, maximal unitarity is the natural continuation of the direct extraction procedures for one-loop integrals. Maximal unitarity relies on a unitarity comptatible basis of linearly independent integrals and the fact that basis elements necessarily involve nontrivial tensorial numerators that contaminate the cuts is a difficulty. However, steps towards developing a general basis of planar two-loop integrals were recently taken in [44,45]. The idea of maximal unitarity is to isolate individual integrals in a basis by localizing their integrands onto global poles and thereby cut as many propagators as possible. The loop amplitude in question falls apart and becomes a product of its constituent trees which upon integration along complex contours yield very compact final expressions whose simplicity is diametrically opposite to that suggested by the Feynman diagram approach. Until now, integral bases for four-particle amplitude contributions have been constructed by reduction to a set of master integrals via integration-by-parts identities. Maximal unitarity has been successfully applied in hepta-cuts of planar and non-planar double box amplitude contributions with and without massive external legs in generic gauge theories [46][47][48][49][50][51][52]. Previously octa-cuts and hepta-cuts of two-loop amplitudes were addressed in maximally supersymmetric N = 4 super Yang-Mils theory [42,43].
The current status of maximal unitarity is that it has been demonstrated to work for univariate residues. Indeed, hepta-cuts of two-loop amplitudes generate algebraic curves [48,63] and the unfrozen degree of freedom on the maximal cut is reflected by a left-over one-dimensional contour integral that computes residues. It is common that amplitudes at higher loops or below the leading singularity at the two-loop level require several complex variables to parametrize a solution of a unitarity cut. This necessarily leads to residues of higher-dimensional differential forms. It is of obvious theoretical and practical interest to extend the univariate unitarity method to general amplitude contributions. The aim of this paper is therefore to develop a systematic way of determining master integral coefficients in an integral basis using multivariate residues in generalized unitarity cuts that define multidimensional algebraic varieties.
The paper is organized in three main parts in the following way. In section 2 we introduce the mathematical prerequisites of multivariate residues and computational algebraic geometry. Then, in section 3 we apply the formalism to the planar triple box and derive a general formula for the master integral coefficients. Finally, in section 4 we obtain explicit results for the three-loop triple box in any renormalizable gauge theory for independent helicity configurations of four external gluons. We choose this primitive amplitude because it was calculated short time ago, yet it is a highly nontrivial test of our method.

Multivariate Residues
Motivated by the discussion in the introduction we briefly review residues for several complex variables. We focus on the computation of degenerate multivariate residues and the global residue theorem. For mathematical concepts, the references are the textbooks by Hartshorne [79], and by Griffiths and Harris [80].
We are working with n complex variables, namely, z ≡ (z 1 , . . . , z n ). First we study the local properties of a residue. Consider a residue at (z 1 , . . . , z n ) = (ξ 1 , . . . , ξ n ) ≡ ξ. Let U be the ball ||z − ξ|| < ǫ and assume that the functions f 1 (z), . . . , f n (z) are holomorphic in a neighborhood of the closureŪ of U , and have only one isolated common zero, ξ in U . Let h(z) be a holomorphic function in a neighborhood ofŪ . Then for the differential form, the residue at ξ regarding the function list {f 1 , . . . , f n } is defined to be, where the contour Γ is defined by the real n-cycle Γ = {z : |f i (z)| = ǫ i } with the orientation specified by the differential form d(arg f 1 ) ∧ · · · ∧ d(arg f n ). The value of a local residue does not only depend on the location, but also the ordering in the function list {f 1 , . . . , f n }. However, to simplify the notation, we may rewrite the left hand side of eqn. (2.2) as Res ξ (ω). The value of a residue is invariant under a nonsingular change of complex coordinates. And we can rescale the function list, f i (z) → α i (z)f i (z), where α i (z)'s are holomorphic functions and α i (ξ) = 0. Although the contour is changed, Res {f 1 ,...,fn},ξ (ω) = Res {α 1 f 1 ,...,αnfn},ξ (ω), from the fact that ω is closed and Stokes' theorem.
It is easy to prove that if locally h(z) is generated by {f 1 , . . . , f n }, i.e., h(z) = a 1 (z)f 1 (z) + · · · + a n (z)f n (z) , where the a i (z)'s are holomorphic functions in a neighborhood of ξ, then The fact comes from Stokes' theorem [80]. We define the local residue to be non-degenerate, if the Jacobian of {f 1 , . . . , f n } at ξ is nonzero, In this case, the value of residue is simply [80], There is another type of residues which can be calculated straightforwardly, factorizable residues. The definition is that each f i is a univariate polynomial, namely, f i (z) = f i (z i ). In this case, the contour in eqn. (2.2) is factorized to the product of univariate contours, Then, by using the univariate residue formula n times, we get the value of the residue. However, in general, a residue is neither non-degenerate nor factorizable. Usually it is not convenient to use the definition (2.2) to calculate a residue directly. Hence we need new techniques for such computation.

Calculation of (degenerate) local residues
We restrict our computation to cases in which f 1 , . . . , f n are polynomials in (z 1 , . . . , z n ). This is sufficient for all unitarity calculations in quantum field theory. When all f i 's are polynomials, we can use the powerful tool of computational algebraic geometry to obtain the residue.
The key is to transform the polynomial list {f 1 , . . . , f n } to a new polynomial list {g 1 , . . . , g n } such that the new residue becomes factorable. So we recall the transformation law, Theorem 1 (Transformation law). Let I = f 1 , . . . , f n be the zero-dimensional ideal generated by {f 1 , . . . , f n } and J = g 1 , . . . , g n be a zero-dimensional ideal such that J ⊂ I. So g i = a ij f j , where the a ij 's are polynomials. Let A be the matrix of a ij 's, then for residues at ξ, Res {f 1 ,...,fn},ξ h(z)dz 1 ∧ · · · ∧ dz n f 1 (z) · · · f n (z) = Res {g 1 ,...,gn},ξ h(z)dz 1 ∧ · · · ∧ dz n g 1 (z) · · · g n (z) det A . (2.8) For the proof of this theorem, we refer to [80]. This theorem holds for both nondegenerate and degenerate residues.
In practice, we use Gröbner basis method to find a list of {g 1 , . . . , g n } such that g i (z) = g i (z i ). To obtain g i , we calculate the Gröbner basis of I = f 1 , . . . , f n in Lexicographic order with the ordering z i+1 ≻ z i+2 ≻ · · · z n ≻ z 1 ≻ z 2 ≻ · · · ≻ z i . Since this lexicographic order would eliminate variables z 1 , . . . , z i−1 , z i+1 , . . . , z n , there would be a polynomial in the basis which only depends on z i . Then we call this polynomial g i . Repeating this process n times, the desired {g 1 , . . . , g n } is obtained. After applying the transformation law, the residue becomes factorable so it can be calculated.
We demonstrate this algorithm by several examples, There is a degenerate residue at (0, 0). Note that three factors z 1 , z 1 + z 2 and z 1 − z 2 vanish at the residue so the Jacobian is zero.
For our residue calculation of the three-loop triple box, we frequently deal with the degenerate residues demonstrated in Example 1 and 2, as (1) more than 2 factors vanish at the same point (2) one of the vanishing factor has the power larger than 1. These degenerate residues are all calculated by Gröbner basis method. The computation is automated by a program 1 powered by the algebraic geometry software Macaulay2 [78].
Furthermore, when n = 2, if three different factors f 1 , f 2 and f 3 vanish at the same point ξ, we have three different residues at ξ, namely, (2.13) In general, the values of three residues are not the same. However, it is clear that for a large number of classes of f i 's, the sum of three residues is zero.
Similarly, it is easy to check that the identity holds when one of the f i is a power of a linear function or an irreducible quadratic polynomial. In our triple box residue computation, we will use this lemma and its generalized version to pick up independent residues.

Global Residue Theorem
It is well known that for the univariate complex analysis, the sum of residues of meromorphic differential form on CP 1 is zero, by the residue theorem. We review the multivariate version, the global residue theorem, in this subsection. This theorem relates the residues in C n and residues at infinity, so it is crucial for the study of the global structure of residues.
Theorem 2 (Global residue residue, GRT). Let M be a n-dimensional compact complex manifold, D 1 , . . . , D n be divisors of M . Assume that the intersection S = D 1 ∩ D 2 · · · ∩ D n only consists of discrete points. ω is a holomorphic differential n-form on M −D 1 ∪D 2 · · ·∪ D n . Then regarding the divisor list {D 1 , D 2 , . . . , D n }, P ∈S Res P (ω) = 0 . (2.18) A proof for this theorem can be found in [80]. Here are some explanations for this theorem: locally, the complex manifold M is C n , so we can set up a local coordinate system (z 1 , . . . , z n ). Also, locally each divisor D i is an analytic hypersurface defined by f i (z 1 , . . . , z n ) = 0. Then the residue regarding {D 1 , . . . , D n } at P is defined to be the local residue at P regarding {f 1 , . . . , f n }. From previous discussion, the value of the residue is independent of the coordinate choice or local rescaling of the f i 's. So it is well-defined.
Although we may start our analysis from C n , the manifold C n is not compact so GRT does not apply directly. So we consider the projective space CP n with the homogenous coordinates [w 0 , w 1 , . . . , w n ]. CP n is covered by the open sets U i , where U i is defined to be the set {[w 0 , w 1 , . . . , w n ]|w i = 0}. C n is embedded inside CP n as U 0 , The points with the coordinate w 0 = 0 are called points at infinity, which form the space CP n−1 . CP n is compact and on which we can use GRT.
Example 3. Consider n = 2 and the differential form We study this form on CP 2 . On the patch To use GRT, we may define the two divisors, However, it is also possible to define two divisors D ′ , [0, 0, 1]}, so again there are two residues. Again GRT holds, So in summary, for ω, apparently there is one residue at the origin, two residues at infinity. By using GRT twice for different choices of divisors, we can see that values of the two residues at infinity are both the opposite to that at origin.
This simple example demonstrates GRT for multivariate residues. In contrast to the univariate case, there are many different residue identities for one differential form, by choosing different divisors. In this paper, we will use GRT to study the residues at infinity.

Maximal Cuts of Planar Triple Box Integrals
We consider color-stripped four-gluon scattering amplitudes at three-loops in gauge theories with SU(N c ) symmetry group. The complete amplitude receives contributions of a large number of integral topologies, for instance the three-loop ladder with one or two crossed boxes and the planar tennis court. Our focus is here on a primitive amplitude of modest complexity, the planar triple box, although our method should generalize to other examples in a straightforward manner.
The dimensionally regularized Feynman scalar integral for the four-point planar triple box with massless kinematics on both internal and external lines reads where the denominators or inverse propagators according to fig. 1 are given by All external momenta are outgoing and summed as K i 1 ···in = k i 1 + · · · + k in . Our definition of the triple box integral is equivalent to that in [55] up to a linear transformation of each of the three loop momenta, In general, (3.1) has a nontrivial polynomial numerator function P(ℓ 1 , ℓ 2 , ℓ 3 ) inserted and is in that situation referred to as a tensor integral, although all Lorentz indices are properly contracted. We specialize to four dimensions and therefore only reconstruct the master integral coefficients to leading order in the dimensional regulator ǫ, leaving the general case to future work. The ten inverse propagators generate a polynomial ideal I = D 1 , . . . , D 10 and the deca-cut equations, i.e. all internal propagators are placed on their mass shell, define a two-dimensional algebraic variety (or four-dimensional real surface) which is the zero locus of I, usually denoted Z(I). In our notation, the solution set is By primary decomposition of the ideal it can be proven that the variety is reducible and the algebraic set can be decomposed uniquely into a union of fourteen components which are in one-to-one correspondence with the inequivalent deca-cut solutions [55,56], The branches come in complex conjugate pairs, although, technically speaking, identification by complex conjugation presumes reality of the external momenta. We label the solutions S 1 , . . . , S 14 such that S i and S i+7 form a pair for i = 1, . . . , 7. Note that S i is not a Riemann sphere because of the dimension. This number of solutions is actually expected from general considerations of the Lloop ladder topology with chiral vertices. Indeed, such a diagram splits into 2(L + 1) three-point trees on the maximal cut with (3L + 1) propagators placed on-shell. Let N L be the number of solutions at L loops with N 1 = 2 and N 2 = 6. We can consider the distribution of holomorphically and anti-holomorphically collinear vertices and easily realize that N L = 2(2 L − 1). In particular, N 3 = 14, N 4 = 30, N 5 = 62 et cetera.
We proceed from here in three steps: first we solve the cut equations, then we localize the scalar master integral onto the on-shell branches to expose the composite leading singularities of the loop integrand, and finally we derive unique projectors for all master integrals that are consistent with integral reduction identities from parity-odd terms and total derivatives which integrate to zero.

Parametrization of On-Shell Solutions
Uniformity of the dimension across all deca-cut branches implies that we can parametrize the on-shell surfaces with two variables (z 1 , z 2 ) ∈ C 2 . We exploit experience from previous calculations at one and two loops and choose convenient normalizations and basis elements such that the global poles of the integrand are directly exposed and as many of the decacut constraints as possible are linearized. Our parametrization of the three independent loop momenta with two-component Weyl spinors follows that of [55] and reads All parameters are complex valued. In advance of calculations below we define a frequently used ratio of Mandelstam invariants, In order to compensate for the change of variables from momenta to parameters in the triple box integral we in principle have to include the three Jacobian determinants .
They are however constant as an immediate consequence of linearity of the parametrization of the loop momenta and can thus be disregarded in the augmented deca-cut in what proceeds. Let us write out the deca-cut equations using the parametrizations of ℓ 1 , ℓ 2 and ℓ 3 . It is a straightforward task and we immediately derive all equations involving only one loop momentum and one extenral leg, Moreover, we have For generic external kinematics, we readily get the trivial constraints α 1 = α 2 = 0 and β 1 = β 2 = 0 together with α 3 α 4 = 0 and β 3 β 4 = 0. It then remains to consider the inverse propagators on the two middle rungs of the triple box. We omit for brevity the full expansions in parameter space, and display only equations simplified on the constraints for ℓ 1 and ℓ 2 above, The complete set of local solutions to these equations were reported in [55]. We quote the result for S 1 , . . . , S 7 below for completeness. Also, the complex conjugates denoted with primes can be constructed easily using the following relations valid for each pair, On the maximal cut, all eight vertices have three massless legs attached. We point out that the fourteen on-shell branches are in one-to-one correspondence with the valid kinematical configurations of three-point trees in the maximally cut planar triple box such that no external legs are neither holomorphically nor antiholomorphically collinear for generic momenta. All diagrams not related by complex conjugation are included in Appendix B. Table 1. Values of the coefficients in the outermost loop momenta ℓ 1 and ℓ 2 written in terms of the two unfrozen parameters (z 1 , z 2 ) ∈ C 2 . Table 2. Solutions for loop momentum ℓ 3 .

Composite Leading Singularities
Loop momenta that solve simultaneous on-shell contraints are complex valued for general external kinematics. As a direct conseqeuence, the traditional unitarity cut procedure involving several Dirac delta functions that enforce on-shell conditions would yield a trivial equation reflecting that their support is zero and this is not sufficient to extract the integral coefficients. This leads to an obstacle already for quadruple cuts at one loop in the work by Britto, Cachazo and Feng [8], later clarified by Kosower and Larsen [46]. Instead of a integration region that comprises real Minkowski space, loop integrals should really be reinterpreted in complex momenta, ℓ µ i ∈ C 4 , evaluated along the real slices Im ℓ µ i = 0. Generalized unitarity cuts that are also appropriate to complex solutions may then be realized by replacement of contours by a higher-dimensional surface that is topologically equivalent to a torus encircling a global poles of the loop integrand. Recall the transformation of localization property of multidimensional contour integrals, which we immediately recognize as the natural generalization of the Dirac delta function to several complex variables once we for ξ ∈ C n apply the substitution . (3.12) Notice that, crucially, analyticity of the integrand is maintained in virtue of absence of absolute values in the Jacobian. As a consequence we are able to carry out successive contour integations. Let us apply this strategy to compute the localization of the triple box scalar integral onto the fourteen on-shell branches and begin to explore the composite leading singularities of the integrand. Each integral must be treated individually. However, due to the appearance of complex conjugate pairs of solutions to the deca-cut equations, we expect at most seven distinct remaining two-dimensional complex integrals. The general object of our interest is where T α,i × T β,i × T γ,i is the twelwe-dimensional leading singularity cycle associated with solution S i and the prefactors are given in (3.8). This multivariate residue is nondegenerate, hence performing ten of the integrals by means of the localization property for multidimensional contour integrals (2.6) leaves us with a two-dimensional contour integral and the deca-cut Jacobian J i defined via (3.14) To set the stage we derive the expression for J 1 in full detail. In this case we keep integrals over α 3 and β 3 . Other choices yield equivalent end results as long as the coordinate transformation to z 1 and z 2 is nonsingular. Indeed, this change of variables introduces of course yet another Jacobian which is trivial to compute. Anyway, the Jacobian that arise upon cutting the ten propagators on-shell is block diagonal and decomposes into three determinants. In particular, we denote the Jacobian from the six deca-cut equations involving only either loop momentum ℓ 1 or ℓ 2 by J A , in a notation that should be self-explanatory. It is then easy to get The next step is to integrate over the four constraints involving ℓ 3 and evaluate the determinant which we call J B . Prior to insertion of the explicit cut solution, this expression expands to a slightly complicated form which we do not include here, Finally we apply the coordinate transformation (α 3 , β 3 ) → (z 1 , z 2 ) and obtain the integral .
We repeated the calculation for the remaining thirteen deca-cut solutions and verified that complex conjugate pairs indeed have the same Jacobian up to an ambiguous overall sign due to antisymmetry of determinants, i.e.
for i = 1, . . . , 7. In summary, we found the following relatively simple integrals 20) , (3.21) , (3.23) , (3.24) We emphasize that the number of polynomials in the denominators and hence residues to consider differs from branch to branch, which is not the case in the two-loop double box, reflecting a more complicated global topological picture of the underlying algebraic variatey. In general, these contour integrals have a nontrivial polynomial numerator which is evaluated on the particular branch. Although the Jacobian computations are not affected by the presence of a numerator function, the analytic structure of the integrands may change significantly. Notice that the cut integrals in solution S 7 and its complex conjugate are particularly simple because they separate onto univariate contours with factorizable residues which can be computed iteratively as in (2.7).

Global Poles and Augmentation of Residues
We expand the planar triple box contribution to the four-gluon three-loop amplitude onto an integral basis of three master integrals, say, where the ellipses cover terms with less than ten propagators which cannot be probed by maximal unitarity. It of course possible to construct other equivalent integral bases for the ladder topology via integration-by-parts identities, however this choice allows to directly compare with the results of Badger, Frellesvig and Zhang. The analytic expression of the scalar master integral I P 10 [1] is known in dimensional regularization [74]. Our problem is now reduced to extract the three master integral coefficients to O(ǫ 0 ) from the augmented deca-cut of (3.26). It is customary to strip off universal prefactors and define dimensionless coefficientsĈ i by where the four-gluon tree amplitude is given by the Parke-Taylor formula, for instance for helicities − − ++, Notice that we form now on leave helicity labels on coefficients implicit. Then we can write down the augmented deca-cut of (3.26) schematically. The triple box amplitude contribution falls apart into a product of eight on-shell three-particle amplitudes such that where the tree-level data evaluated in the deca-cut kinematics defines fourteen multivariate Laurent polynomials, All coefficients with |i| or |j| > 4 vanish in all gauge theories even though higher order terms are allowed in renormalization [55]. It follows from the explicit solutions that the Laurent series for S 1 and its complex conjugate S 8 contain terms with negative powers of z 1 and z 2 in the numerators. Furthermore, parameters in S 2 , . . . , S 4 and their complex conjugates contain negative powers of z 1 . These will also give rise to global poles. We therefore accordingly construct 14 residue loci from all combinations of two vanishing polynomial factors in the deca-cut Jacobians and the possible extra denominators z 1 , z 2 . Here we list all residues for finite values of (z 1 , z 2 ). The residues are labeled by the number of branch in Roman numerals, and its order in this branch. Multivariate residues also depend on the choice and order of the denominators. So in (...) we give the location of the residue, and in {. . .} we list the denominators vanishing on the residue. There are also residues at infinity. However, we prove that for the maximal cut of the triple box, all residues at infinity are related to poles at finite loci by the Global residue theorem (GRT). So we do not list residues at infinity and we give the proof in Appendix A.
Second, the 14 branches intersect with each other and many residues are corresponding to the same value of the loop momenta. For example, it is clear that, at the residues I1 and II5, the ISP values are the same, (3.31) So they correspond to the same point in loop momenta space. Then we can remove residue II5. Sometimes, the parametrization apparently blows up at one residue. Then we can consider the projective space of loop momenta and study the coincident residues. For example, residues I5 and IV3 correspond to the same point in the projective loop momenta space. Exhausting all these obvious coincident residues, we further remove 35 residues.
Third, there are more subtle coincident residues. For example, on the branch S 1 , the limit z 1 → −1, z 2 → 0 does not correspond to one point in the projective loop momenta space, since the limit depends on the path approaching (−1, 0). Explicitly, we can check that for the small contour defined in (2.2) with the denominators choice of I2, the corresponding loop momenta approach that of residue V1. Similarly, (I3, V3), (VIII2, II2), (VIII3, XII3) are also coincident residues in this pattern. So we can further remove 4 residues.
Finally, we have 23 independent residues and choose them to be, Then we calculate the values of all terms in the numerator on these residues. The nondegenerate residues can be calculated directly by the Jacobian matrix, while the degenerate residues are calculated by the transformation law (2.8). We use a program powered by the algebraic geometry software Macaulay2 [78]. Actually we calculated the values on all 64 residues, and explicitly verified the relation (2.14) and that for each coincident residue pair, the values are the same. Now that we have a minimal basis of triple box residues we can transform (3.29) into an algeberaic equation for the master integral coefficients. We simply expand the master integrals, namely I   where the values are listed in the same order as (3.32).

Integral Reduction Identities
Consistency of the unitarity procedure of course requires that vanishing Feynman integrals must also have vanishing deca-cuts. However, promotion of real slice integrals to arbitrary multidimensional contour integrals implies that certain integral relations eventually fail to hold. Therefore it is necessary to constrain the integration contours and thus demand that any valid integral identity is preserved under replacement of contours, which is to say, I 1 = I 2 =⇒ cut(I 1 ) = cut(I 2 ) . (3.36) The origin of the integral reduction identities used for projection onto master integrals is simple. In order to understand this better let us briefly consider the general integrand numerator polynomial which can be parametrized in terms of seven irreducible scalar products and naturally splits into spurious and nonspurious parts, where the additional orthogonal direction is represented by The maximum powers of the coefficients are restricted by renormalizability requirements and completion of the integrand reduction is achieved by multivariate polynomial division with respect to a Gröbner basis constructed form the ten inverse propagators using the Mathematica package BasisDet [56] written by one of the present authors. The final integrand then consists of 199 spurious and 199 nonspurious terms. Although spurious terms are of course not present in integrated expressions, they play a vital role at the level of the integrand and for the winding numbers. We require that all parity odd vanishing integral identities continue to hold after pushing integration contours into (C 4 ) ⊗3 . This leads to nontrivial constraints on the winding numbers. Even powers of spurious ISPs are reducible and can be expressed in terms of other scalar products by means of Gram determinant identities for four-dimensional momenta. Another equivalent strategy is to identify the full variety of Levi-Civita insertions that appear in integral reduction, after using momentum conservation.
The fact that total derivatives vanish upon integration allows us to add such terms to the integrand and thereby produce a vast set of relations among integrals known as integration-by-parts identities of the form In order to ensure validity of the unitarity method each nontrivial identity gives rise to a constraint requiring that the procedure yields a vanishing coefficient for the additional term as this is not true automatically contour by contour. In the four-point massless planar triple box this amounts 199 integral identities corresponding to the nonspurious part of the integrand. We include some of those identities for illustration purposes in Appendix C. For most families of Feynman integrals the relations are often quite involved to obtain. We refer to the public computer codes FIRE [76] and Reduze [77], for more information on how to generate the identities in practice.

Unique Master Integral Projectors
Based on the discussion of the previous sections, we compute the residues from the deca-cut to resolve all constraints from integration-by-parts identities and vanishing of all spurious terms upon integration, and organize them as a homogeneous system of equations with a matrixM of size 398 × 23. It turns out that rankM = 20. In other words, the integration contours are subject to only 20 constraints in order to yield a valid unitarity procedure.
Overall we are therefore left with three independent contour weights that are not fixed by integral reduction consistency requirements, exactly matching the number of planar triple box master integrals. This pleasant freedom ensures that we can derive projectors or master contours that normalize one master integral to unity and set the remaining two to zero. We therefore extend the matrixM to a 401 × 23 matrix M with the residues of the master integrals and write down equations for the projectors, The associated master contours M i can be constructed explicitly as linear combinations of weighted infinitesimal toroidal surfaces encircling the 23 global poles. According to (3.32), each master contour only receives contributions from a small subset Λ of the 14 on-shell branches after removal of redundant residues. Therefore we can apply the decomposition M i = k∈Λ M i;k for Λ = {I,II,III,V,VIII,IX,XI}. In this notation, weights are kept implicit in the contours. The master integral coefficients can then be written schematically in the very compact form This formula completes our derivation of the triple box master integral coefficients.

Examples
In this section we apply the master integral coefficient formulae to the triple box contribution to three-loop four-point gluon amplitudes with specific helicity configurations and massless kinematics. We only consider deca-cuts in the s-channel as contributions in the t-channel can be obtained in a completely similar manner by accounting for cyclic permutation of external legs. The starting point for the computation is the intermediate state sum over the product of eight tree-level amplitudes that arise by cutting the numerator function N on-shell on each branch, The expression is summed over all possible internal helicity states and distinct configurations of gluons, fermions and scalars inside the loops. It may be obtained quite easily in N = 0, 1, 2, 4 supersymmetric Yang-Mills theories by superspace techniques [65,66] or simply by direct computation in a generic theory from the distributions of holomorphic and antiholomorphic vertices (see Appendix B) and all possible flavour configurations. For now we exploit that the triple box tree-level data has already been calculated for n f fermion and n s complex scalar flavours in the adjoint representation and refer to [55] for more details. We examine the alternating helicity configuration − + −+. In this case, the 23 independent residues computed from the Laurent expansions of the products of tree-level amplitudes explicitly read where the integrand is understood to be evaluated on the branch on which the residue R i actually resides. The new variables r 1 and r 2 introduced in (4.2) are defined in the following way, In a supersymmetric theory we may rewrite n f = N and n s = N − 1 where N counts the number of supercharges. With this in mind, the residues were resolved into a form with coefficients that vanish for N = 4 and also N = 2. Summation of all residues with weights according to the three master contours now yields the normalized master integral coefficients,Ĉ The case of helicities − − ++ simply giveŝ The coefficients displayed here therefore clearly reproduce the known result in N = 4 super Yang-Mills theory [75] and also agree with the general expressions for any number of adjoint fermions and scalars [55]. We finally verified the last independent helicity configurations, namely − + +−.

Conclusion
In this paper we have extended four-dimensional univariate maximal unitarity [46] to cuts that define multidimensional algebraic varieties in any gauge theory. As an example, we studied the application to three-loop amplitudes. In maximal unitarity one cuts the maximum number of propagators on-shell by replacement of real slice integrations by multidimensional contours encircling the global poles of the loop integrand. Each residue comes with a weight or winding number in order for the procedure to conform with integral reduction identities from integration-by-parts relations and parity-odd terms that vanish upon integration.
The technique was demonstrated explicitly for the planar triple box with four external massless legs. We obtained unique and very simple projectors for all master integrals from just 23 finite multivariate residues of the maximally cut loop integrand. We proved by the Global Residue Theorem that all residues at infinity are linear combinations of finite residues and therefore not needed in the computation. Also, we worked out master integral coefficients for all independent helicity configurations in any gauge theory and found exact agreement with recently obtained results [55].
Our method is completely general and we expect it to apply to multivariate residues in more than two complex variables, for instance in the L-loop ladder topology and similar planar and nonplanar diagrams. It also establishes an initial foundation for extracting master integral coefficients below the leading singularity at the two-loop level and therefore provides a step towards full automation. A complete calculation also requires terms of O(ǫ) beyond four-dimensional cuts.
We also give directions for research in the future. It would be very beneficial to better understand the simplicity of the projectors at two and more loops. Guided by recent progress in integrand-level reduction [57,64], we find it very interesting to generalize the maximal unitarity method to D = 4−2ǫ dimensions where amplitudes are cut into products of six-dimensional trees. In that case, the polynomial ideal generated by a set of propagators is prime and hence there is only one branch of the unitarity cut. This will maybe make it easier to resolve the redundancy among residues. Scattering processes with more than four external particles are also important to examine. Higher-point integrals suffer from a more complicated set of reduction identities due to dependence on a broader variety of kinematical invariants. Planar and nonplanar double box integrals with five massless external legs probably hint the simplest extension of maximal unitarity in this direction because their hepta-cut equations are similar to those of the four-point case. As a general mathematical tool we believe that multivariate residue calculations will also prove useful broadly in the the study of supersymmetric Yang-Mills amplitudes and evaluation of master integrals. We also expect that computational algebraic geometry can offer important insight in the topological information of algebraic varieties from unitarity cuts such as degeneracies under specific kinematics. Hopefully some of these questions will be addressed soon.

A Residues at Infinity
In this Appendix, we use the global residue theorem (GRT) (2.18) to prove that for the planar triple box, all residues at infinity are linearly related to residues at finite loci.
Since the branches S 8 , . . . , S 14 are complex conjugate of S 1 , . . . , S 7 , we only present the residues at infinity for S 1 , . . . , S 7 . We start from the simplest cases S 5 , S 6 and S 7 , then S 1 and S 2 , and finally the complicated cases S 3 and S 4 .
We extend the space C 2 : (z 1 , z 2 ) to CP 2 : [w 0 , w 1 , w 2 ] as z 1 = w 1 /w 0 and z 2 = w 2 /w 0 . The points in CP 2 with w 0 = 0 are called points at infinity. • Branch S 5 . The denominators in the Jacobian are z 1 , 1 + z 1 − z 2 , z 2 . Furthermore, the differential form may blow up at w = 0. So we list the vanishing polynomials in projective variables, It is clear that w 1 and w 0 both vanish at [0, 0, 1], w 0 + w 1 − w 2 and w 2 vanish at [0, 1, 1], w 2 and w 0 vanish [0, 1, 0]. So there are 3 residues at infinity. Each of these residues are linear combinations of the residues at finite loci, by GRT. For example, for the residue with the denominators, we can choose two divisors on CP 2 , We find it is helpful to illustrate the residues at infinity by a sketch diagram in fig. 2. (e). We use lines (curves) to represent vanishing polynomials and intersection points (black dots) for residues at infinity. Residues at finite loci are not shown. We call such a diagram, infinity diagram. For branch S 5 , we use a horizontal line for the polynomial w 0 , and three vertical lines for w 1 , w 0 + w 1 − w 2 and w 2 . There are three intersecting points, which are three residues at infinity. To use GRT for [0, 0, 1], we split the diagram to two components, one with w 1 , the other with w 0 w 0 +w 1 −w 2 and w 2 . These two components correspond to two divisors. Because the two components have no other intersecting point at infinity, by GRT, the residue [0, 0, 1] should be a linear combination of residues at finite loci. This diagram makes the proof clear.
• Branch S 6 . The vanishing polynomials in projective variables are, and there are three residues at infinity. Its infinity diagram, fig. 2. (f), has the same structure as that of S 5 , and by the same proof, all residues at infinity are linear combinations of residues at finite loci.
• Branch S 7 . The vanishing polynomials in projective variables are, The first two polynomials and w 0 vanish at [0, 0, 1], while the third and fourth and w 0 vanish at [0, 1, 0]. There are two residues at infinity and both are degenerate, which can give four values depending on the choice of denominators. The infinity diagram is shown in fig. 2 (g). Again, the topology of this diagram, is a tree. So for any residue at infinity, no matter how to choose the denominators, we can always split the diagram into two and the two denominators are in different components. Then GRT shown that all residues at infinity can be ignored.
• Branch S 1 . The vanishing polynomials in projective variables are, The intersecting structure is more complicated. However, the infinity diagram, shown in fig. 2. (a), is still a tree, so by using GRT several times, we can shown that all residues at infinity can be ignored.
• Branch S 2 . The vanishing polynomials in projective variables are, Similarly, the infinity diagram, shown in fig. 2. (b), is a tree, so all residues at infinity can be ignored.
• Branch S 3 . The vanishing polynomials in projective variables are, w 1 , w 0 + w 1 , w 2 , w 2 0 + w 0 w 1 + χw 0 w 1 − χw 1 w 2 , w 0 . (A.7) The structure of residues at infinity of this diagram is complicated. Note that the fourth polynomial is quadratic, and it has two intersection points with the hypersurface {w 0 = 0}. The first, second, fourth and fifth polynomials vanish at [0, 0, 1]. On the other hand, the third, fourth and fifth polynomials vanish at [0, 1, 0]. So the corresponding infinity diagram, shown in fig. 2. (c), has the topology of a loop. Hence the direct proof by GRT does not work for this branch, and residues at infinity are not linear combination of residues at finite loci on this branch.
However, we show that the residues at infinity are actually the sum of residues on S 3 and also other branches. It is clearly that by GRT, if we know all the values of the residue at [0, 1, 0], we get all the values of residues at [0, 0, 1]. At [0, 1, 0], there are three residues, corresponding to, {w 0 w 2 , w 2 0 + w 0 w 1 + χw 0 w 1 − χw 1 w 2 } {(w 2 0 + w 0 w 1 + χw 0 w 1 − χw 1 w 2 )w 0 , w 2 } and {w 2 (w 2 0 + w 0 w 1 + χw 0 w 1 − χw 1 w 2 ), w 0 }. By Lemma 1 (2.14) (generalized version), the three residues sum to zero, so we only need the first two. Furthermore, for the second one, the two polynomials for the loop in the infinity diagram are combined together. So it is possible to split the diagram into two components intersecting at only one point. Then this residue is a combination of residues at finite loci. Finally we only need to consider the first one. Explicitly, we find that this residue and the residue II4 (on branch S 2 ) correspond to the same point in the projective loop momenta space. The values of this residue for all integrand terms match these of residue II4. So this residue at infinity is not needed. Therefore, no residue at infinity is needed for this branch.
• Branch S 4 . The vanishing polynomials in projective variables are, w 1 , w 0 + w 1 , w 2 , w 2 0 + w 0 w 1 + χw 0 w 1 + w 1 w 2 + χw 1 w 2 , w 0 . (A.8) Again, the infinity diagram, shown in fig. 2. (d), is a loop and so the direct proof by GRT does not work. However, by the same analysis as that for S 3 , the residues at infinity are actually the sum of residues on S 4 and also other branches. No residue at infinity is needed for this branch.

B Kinematical Configurations of the Planar Triple Box
We include here valid kinematical configurations of three-point trees in the maximally cut planar triple box such that no external legs are neither holomorphically nor antiholomorphically collinear for generic momenta. The diagrams are in one-to-one correspondence with the on-shell solutions and follow their enumeration described in the main text. For brevity we only depict contributions that are not related to each other by parity conjugation. In our conventions, black and white blobs denote MHV and googly-MHV vertices respectively.

C Planar Triple Box Integration-By-Parts Identities
We provide below a small subset of the four-dimensional integration-by-parts identities used for the reduction onto master integrals of renormalizable Feynman integrals with planar triple box topology. The relations are for simplicity stated in terms of irreducible scalar products (ISPs) defined by x ij ≡l i · v j with momentum basis v = (k 1 , k 2 , k 4 , ω) so that the master integrals can be written I (3 + χ(3 + χ))s 5 12 I P 10 [x 13 ] + · · · (C.8) Ellipses denote truncation at the maximal number of propagators. We also have to consider integral reduction identites for terms in the integrand that contain products of spurious scalar products x 14 and x 24 . A few of these relations read The remaining relations, which we do not include here, display the same simplicity. Many high-rank integrals are actually reducible below ten propagators, which means clearly means that they should vanish on the maximal cut.