Massive Nonplanar Two-Loop Maximal Unitarity

We explore maximal unitarity for nonplanar two-loop integrals with up to four massive external legs. In this framework, the amplitude is reduced to a basis of master integrals whose coefficients are extracted from maximal cuts. The hepta-cut of the nonplanar double box defines a nodal algebraic curve associated with a multiply pinched genus-3 Riemann surface. All possible configurations of external masses are covered by two distinct topological pictures in which the curve decomposes into either six or eight Riemann spheres. The procedure relies on consistency equations based on vanishing of integrals of total derivatives and Levi-Civita contractions. Our analysis indicates that these constraints are governed by the global structure of the maximal cut. Lastly, we present an algorithm for computing generalized cuts of massive integrals with higher powers of propagators based on the Bezoutian matrix method.


Introduction
Contemporary experimental high energy physics is concentrated on the Large Hadron Collider (LHC) at CERN. Our ability to utilize the huge amount of data delivered by the experiment towards further scientific progress relies on a quantitative understanding of all relevant scattering processes in the Standard Model. Otherwise, we are unable to extract signals of new physics from the background. Precise theoretical predictions in Quantum Chromodynamics (QCD) at the LHC require not only amplitudes at leading order (LO) and next-to-leading order (NLO), but also next-to-next-to leading order (NNLO) corrections to comply with the level of accuracy of the data. For some processes, two-loop amplitudes are important already at NLO because the LO terms begin at one loop. The text book approach to perturbative scattering amplitudes is through Feynman rules and diagrams. Although it tracks interactions of particles very intuively and in principle always works, this method suffers from severe computational problems with increasing loop level or number of external legs. The main reason is that the gauge redundancy of the theory introduces virtual intermediate states that are off-shell. Not even very powerful computers are able to deal with many of the phenomenologically interesting processes without new clever ways to attack the problem.
In the last two decades, many attempts to surmount the computational bottleneck have been reported. The lesson is to exploit analyticity and unitarity of the scattering matrix. Analyticity allows for amplitudes to be reconstructed from their singularity structure, whereas by unitarity, residues at the poles factorize onto products of simple amplitudes. Two of the most successful advances are the original unitarity method for loop amplitudes developed by Bern, Dixon and Kosower [3,4] and the Britto-Cachazo-Feng-Witting (BCFW) recursion relations [1,2] for trees. In these works, striking and otherwise completely unexpected structure and simplicity are revealed by virtue of retaining only physical on-shell information in a Lagrangianless setting. In a nutshell, all trees may now be constructed recursively and further fused into loops.
The basic idea of the unitarity method (see also later studies, e.g. refs. [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24]) is to reconstruct the amplitude from double cuts that place internal lines in a given channel on their mass-shell and break it into a product of trees. Many individual contributions share such minimal cuts and are therefore hard to separate. Therefore intermediate algebraic steps are typically needed. In that view, the generalized unitarity method [7,18] is more efficient because several propagators are cut simultaneously and thus fewer integrals are isolated. Thanks to the unitarity method, otherwise unfeasible computations of 2 → 2 massless scattering processes in QCD have been carried out.
In the last couple of years, two-loop amplitudes have received substantial attention in the literature. The integrand-level reduction method of Ossola, Papadopoulos and Pittau (OPP) has been extended to multi loops using computational algebraic geometry, and a general way of classifying high-loop unitarity cut solutions is now available [47,[54][55][56][57][58][59][60][61][62][63][64][65]. These techniques were used by Badger, Frellesvig and one of the present authors to calculate the planar part of the all-plus two-loop five-gluon amplitude in QCD [58] and also demonstrated for the planar triple box [56]. In ref. [66] the unitarity method was applied to two-loop diagrams to determine their integral bases.
Working directly at the level of the integral basis, the maximal unitarity formalism initiated by Kosower and Larsen in ref. [46] has emerged as an extension of the quadruple and triple cut at one loop [7,18]. In maximal unitarity one expands the amplitude in a basis of integrals and seeks to isolate the integral coefficients by finding multidimensional complex integration contours that are uniquely associated with each individual master integral. One of the major advantages of maximal unitarity is that one may circumvent the integrand basis which is typically considerably larger than the integral basis. After the reduction onto master integrals is complete, each coefficient is extracted as a linear combination of residues of the product of trees that arise when the diagram falls apart on-shell. The tree-level data is easy to manage using superspace techniques [67,68]. The leading singularity method [42,43] previously addressed hepta-cuts and octa-cuts at two loops in N = 4 super Yang-Mills theory.
Until now, maximal unitarity has remained relatively unexplored in the nonplanar sector [50]. In this paper, we extend the framework to nonplanar double box integrals with up to four external massive legs. Indeed, inspection of the nonplanar part of the integral basis for, e.g., gg → V 1 V 2 computed in ref. [69] shows that in practice the crossed double boxes constitute most of the nonplanar basis. 1 Remarkably, we find that essentially all features of maximal unitarity observed in the planar sector [46][47][48][49] carry over directly to the nonplanar sector. In particular, in determining projectors for the master integrals, we find that the global structure of the maximal cut seems to govern consistency equations from integration-by-parts (IBP) identities and the number of master integrals. Moreover, we show that the constraints are inherited through chiral branchings between distinct classes of hepta-cut solutions.
The paper is organized in the following way. In section 2 we review the maximal unitarity method and the theory of multivariate residues. In sections 3-5 we respectively parametrize the hepta-cut solutions using mutually projecting kinematics, analyze the global structure of the maximal cut and impose consistency relations in order to uniquely fix the projectors for the master integral coefficients for all kinematically inequivalent configurations. Finally, in section 6, we present an enhanced algorithm to compute degenerate multivariate residues from generalized unitarity cuts and apply the technique to massive integrals with doubled propagators.

Maximal Unitarity
The modern version of the unitarity method relies on the existence of a finite basis of linearly independent master integrals {I i } onto which the amplitude in consideration can be expanded, up to additional rational terms, A L-loop n = i ∈ Basis c i I i + rational terms .
(2.1) Therefore, if the basis integrals are known explicitly in dimensional regularization, calculating an amplitude boils down to determining the rational coefficients {c i }. The trick is to apply generalized unitarity cuts to either site of eq. (2.1), and exploit that the cut amplitude factorizes onto simpler quantities. At the one-loop level, the basis consists of boxes, triangles and bubbles with scalar numerators only. Here the computation has already been fully automated, see refs. [7,18,19,25]. For instance, a box coefficient is isolated by a quadruple cut and thereby becomes the product of the trees at the four corners, evaluated in on-shell kinematics. The on-shell internal lines are complex valued for general external momenta. This implies that the cut prescription in terms of Dirac Delta functions necessarily must be reformulated by means of a multidimensional complex contour integral encircling global poles [46,47].
At two loops and beyond, the situation is more intricate, one of the main reasons being that a minimal integral basis is not yet known. Integrals with numerator insertions are in general algebraically irreducible and the reduction to master integrals inevitably involves IBP identities. As a consequence, multiple contributions contaminate the unitarity cuts. Although maximal cuts for four particles at two and three loops have superficial resemblance to the quadruple cut at one loop, it is also challenging to extract the coefficients because the cut does not localize integrals to a point, but rather an algebraic curve [46,47] or generally speaking, an algebraic surface [51,64].
In the last couple of years, maximal unitarity has been applied to the planar double box with up to four external massive legs [46,48,49] in general theories and to the massless nonplanar double box [50]. Recently, the formalism was also extended to amplitude contributions whose maximal cuts define multidimensional algebraic variaties [51], exemplified for the planar triple box at three loops. Along these lines, the unitarity cut prescription has been extended to accomodate loop integrals with doubled or higher powers of propagators [52].

Multivariate Residues
Inspired by the discussion in the introduction, we start by reviewing basic theory of multivariate residues, with emphasis on computation of nondegenerate residues. We also refer the reader to classical text books by Griffiths and Harris [71] and Hartshorne [72].
Let U = {z ∈ C n : z − ξ < ǫ} for ǫ > 0 be a small ball around z = ξ and assume that f and h are holomorphic maps in a neighborhood of the closureŪ of U . For our purposes, it is in fact sufficient to think of each component of f and h as just being multivariate polynomials of certain degrees. Furthermore, suppose that f −1 (0) ∩ U = {ξ}, i.e. the components of f have exactly one simultaneous zero ξ ∈ U . Then for the meromorphic n-form, the associated residue at z = ξ is computed by an integral over a contour that is topologically equivalent to a torus of real dimension n embedded in C n . In detail we have We remark several elementary properties of the residue. The residue is linear in h, but alternating in the f i 's. Moreover, the value of a residue is invariant under nonsingular complex coordinate transformations. It is not hard to prove by Stokes' theorem that if locally h ∈ I f = f 1 , . . . , f n , where I f is the ideal generated by the f i 's, that is, for holomorphic functions a i in a neighborhood of ξ, then we have For the calculation of a residue we distinguish between three classes of increasing difficult: factorizable, nondegenerate and degenerate residues. If each f i defines a univariate polynomial, i.e. f i (z) = f i (z i ), the contour factorizes onto a product of univariate contours such that the residue can be obtained in a manner that trivially resembles the one-dimensional case, If the f i 's are not univariate polynomials and the Jacobian determinant of f 1 , . . . , f n evaluated at z = ξ is nonzero, We immediately recognize the localization property (2.9) as the obvious generalization of the Dirac delta function to several complex variables once we define . (2.10) In particular, this observation allows us to define generalized unitarity cuts of amplitude contributions that only factorize for complex kinematics. In general, a multivariate residue is neither nondegenerate nor factorizable and we then proceed by means of computational algebraic geometry and use the transformation law and Gröbner basis method as we will explain in section 6. Examples of multiloop unitarity cuts that give rise to degenerate multivariate residues include among others the three-loop planar triple box [51] and integrals with doubled or higher powers of propagators [52]. In the two-loop nonplanar double box computation we will mostly encounter nondegenerate residues.

Parametrization of Hepta-Cut Solutions
The dimensionally regularized Feynman scalar integral for the two-loop crossed box amplitude contribution with possibly massive external legs k 1 , . . . , k 6 distributed across all six vertices is where the inverse propagators are Conventions and momentum flow are shown shown in fig. 1. The Feynman iǫ-prescription has been suppressed as it is irrelevant for our purposes. Generally speaking, this integral will have a nontrivial polynomial numerator function denoted Φ(ℓ 1 , ℓ 2 ) and is in that situation referred to as a tensor integral even though all Lorentz indices are properly contracted. The seven inverse propagators {f i } generate a polynomial ideal I = f 1 , . . . , f 7 and the hepta-cut equations define a complex algebraic curve (or two-dimensional real surface) which is the zero locus of I. In our notation, The curve is generally reducible and the algebraic set S can always be decomposed uniquely into a union of a finite number n of irreducible components which are in one-to-one correspondence with the inequivalent hepta-cut solutions, The existence of such a decomposition can be proven by algebraic geometry and primary decomposition of the polynomial ideal, see ref. [57]. The number of solutions within an integral topology depends on the kinematical configuration, in particular the distribution of massive and massless legs.
In the rest of this paper, we examine four-dimensional amplitude contributions with two-loop crossed box topology with k 5 = k 6 = 0, allowing any other configuration of massive and massless external legs. The Mandelstam invariants are throughout this paper defined as so that momentum conversation can be stated as The four-point massless nonplanar double box was previously studied in terms of residues and multidimensional contour integrals by one of the present authors in ref. [50] and by integrand-level reduction by Badger, Frellesvig and the other author in ref. [54].

Mutually Projected Kinematics
Scattering amplitudes of massless particles are naturally encoded in the spinor helicity formalism by Lorentz invariant innner products of commuting spinors λ α i andλα i . For a momentum k i with k 2 i = 0 we have the representation k αα i = λ α iλα i . It is then possible to define antisymmetric holomorphic and antiholomorphic brackets, and make contact to the Mandelstam invariants, We treat the massive hepta-cut equations for the nonplanar double box using mutually projected kinematics [15,18]. Thereby the spinor helicity formalism and massive momenta become compatible. Given a pair of massive momenta (k i , k j ), the idea is to obtain massless momenta (k ♭ i , k ♭ j ) each of which is the massless projection of one of the massive legs in the direction of the other masslessly projected leg. Here we will consider four-point kinematics with mutually projecting pairs (k 1 , k 2 ) and (k 3 , k 4 ). Other choices are of course also possible. Within each pair we define so that (k ♭,µ j,1 , k ♭,µ j,2 ) are massless momenta by construction. It is easy to verify that (3.10) We streamline notation and define the frequently used quantity which upon identification in eq. (3.9) leads to a quadratic equation whose solutions are It is perhaps more useful from a practical point of view to express the massless projections in terms of the corresponding massive momenta. To that end we let (i,ī) denote a mutually projecting pair and define Then it is straightforward to invert eq. (3.9) with the result On the other hand, the decomposition of the massive legs in terms of a pair of flattened momenta reads It is convenient to introduce short hand notation ρ 12 ≡ ρ 1,1 , ρ 21 ≡ ρ 1,2 and γ 12 ≡ γ j,12 so that ρ 12 = m 2 1 /γ 12 and ρ 21 = m 2 2 /γ 12 , and similarly for the other mutually projecting pair. If m 1 m 2 = 0 we have γ 12 = s 12 and likewise for γ 34 . Our final results can therefore be expressed in terms of the nonzero masses among {m 1 , m 2 , m 3 , m 4 }, γ 12 and γ 34 if respectively m 1 m 2 = 0 and m 3 m 4 = 0, along with two independent Mandelstam invariants, say s 12 and s 14 .
We adopt a loop momentum parametrization of the form, so that ℓ 2 1 = ℓ 2 2 = 0. The various loop spinors are then constructed from the spinors corresponding to the two mutually projecting pairs with general complex coefficients, In this way, massive external momenta are related to holomorphic and antiholomorphic spinors corresponding to their massless projections through 2k µ i = i ♭ |γ µ |i ♭ ]. Expanded explicitly in the basis of four-vectors the two loop momenta read and therefore we are able to eventually fix two of the complex parameters. This freedom amounts sort of a gauge choice. We emphasize that this choice is not necessarily the same for all on-shell branches. In fact, this is not possible in the two-loop crossed box as opposed to the planar case [49].

Four-Mass Hepta-Cut Equations
To begin with we will derive the hepta-cut equations for nonzero external masses in all four corners. Besides the on-shell constraints for ℓ 1 and ℓ 2 which are already satisfied automatically, it thus remains to examine the other five massive hepta-cut equations. Three of them are very simple because they only involve one of the loop momenta. Indeed, it takes little effort to realize that Lorentz products of flattened momenta are needed throughout this calculation. Before we continue let us therefore for completeness derive the necessary expressions. The trick is of course to apply eq. (3.14) and thereby invoke the massive vectors whose contractions are well known. In fact, We also have In addition, we use the same technique to also provide explicit formulae for various contractions of flattened momenta with the external massive legs needed in one of the on-shell equations below, . (3.28) Another rather trivial, but useful, identity is Finally, in what proceeds, we will also encounter the quantities which are complex conjugates of each other for real external momenta as indicated. But actually τ =τ . For completeness, τ can be expanded and reexpressed in terms of the independent kinematical invariants described above in the following way, Let us now return to the hepta-cut equations for loop momentum ℓ 1 . For general masses m 1 = 0 = m 2 we obtain the solution In contrast to the planar double box, there is only one additional nontrivial on-shell constraint for loop momentum ℓ 2 . Rewriting eq. (3.22) in the slightly more suggestive form, The foresight in the choice of parametrization of ℓ 1 and ℓ 2 implies that the coupled on-shell equations are also quite compact actually. It happens that one of them factorizes completely in a symmetric manner, whereas the other can be written  35) or in the slightly more appealing form, .

(3.36)
In the last step we recast the equation by means of the two identities which can be verified through eqs. (3.23)-(3.28) and (3.30) along with expressions forξ 1 andξ 2 given in eq. (3.31), although that task is rather tedious. In order to solve the hepta-cut equations, we decompose the reducible ideal generated by the list of rewritten inverse propagators into an intersection of six prime ideals. Then we compute generating sets that form Gröbner bases over the field of rational functions in each irreducible ideal and obtain the associated zero loci by hand. The six distinct hepta-cut solutions are really three pairs of parity conjugates (S 1 , S 2 ), (S 3 , S 4 ), (S 5 , S 6 ) and each on-shell branch is topologically equivalent to a Riemann sphere, parametrized by a free variable z ∈ C.
We choose to make the behavior under parity conjugation manifest and present the solutions in a symmetric manner such that S 2k−1 and S 2k for k = 1, 2, 3 are related to each other by simply interchanging (ξ 1 , . . . , ξ 4 ) ←→ (ξ ′ 1 , . . . , ξ ′ 4 ). In terms of the parameters , the solutions are The branches are written in terms of independent kinematical constantsξ 1 ,ξ 2 , τ and µ (3.30)-(3.32). For completeness we include the explicit forms of the variables in the four-vector expansions of the two loop momenta for all hepta-cut solutions in appendix A.

Massless Limits
We will also analyze the crossed bouble box with one, two and three massless legs. Integrals for this kinematics are also relevant for higher-multiplicity scattering processes, starting already at five external massless particles. Let us look more closely at the hepta-cut equations and their solutions. We focus on the two momenta in the crossed end of the diagram and assume that m 1 m 2 = 0. The only dependence on m 3 is implicitly through other parameters, e.g. τ and µ. In particular, all on-shell equations and their solutions have the correct limits for m 3 → 0. Moreover, it is clear that µ → 0 for m 4 → 0 so that eq. (3.32) should be replaced by ξ 3 ξ ′ 3 = 0. It can be shown that the number of branches remains six and the explicit solutions follow from the four-mass case once we let µ → 0. The configurations corresponding to this class of kinematics are illustrated in fig. 2.
The situation is slightly more complicated when momenta k 1 and k 2 in the planar end become massless. Here, the massless limits are not smooth and hence they should be The first class of nonplanar double box integrals includes the four-mass case and related massless limits, i.e. three-mass and short-side two-mass with massless legs in the nonplanar end of the diagram. Massless and massive external legs are denoted by single and doubled lines respectively. treated carefully. We see thatξ ′ 1 → 1 as m 2 → 0 and m 1 arbitrary. Also,ξ ′ 1 → 1 + m 2 2 /γ 12 when m 1 → 0. Therefore the equation ξ 1 ξ ′ 1 −ξ ′ 1 = 0 will not give rise to branchings. But if at least one mass among {m 1 , m 2 } is zero,ξ ′ 2 = 0 and we instead get the equation ξ 2 ξ ′ 2 = 0 and eq. (3.31) must be replaced by a pair of solutions, i.e. ξ ′ 2 = 0 and ξ 2 free or vice versa. We solved the hepta-cut equations for this class of kinematics and found eight distinct solutions that can be parametrized as follows,

Classification of Kinematical Solutions
Previous classifications of kinematical solutions in the planar double box with up to four massive legs [47][48][49]53], the nonplanar double box [50] and three-loop triple box [51] with massless external lines suggest that the hepta-cut branches can understood from distributions of vertices. This turns out to be true to some extent in the general four-point two-loop crossed box. Momentum conservation in a three-point vertex forces either square or angle spinor products to align. Phrased slightly differently, the positive or negative chirality spinors of the momenta are collinear, λ a ∝ λ b ∝ λ c orλ a ∝λ b ∝λ c . We choose to depict such vertices as • and • and refer to them as antiholomorphic (MHV) and holomorphic (MHV) respectively. Vertices that involve more than three particles or massive momenta do not have a well-defined chirality.
As a consequence of proportionality of spinors, uninterrupted chains of like-chirality vertices between external legs correspond special kinematics and should therefore be disregarded. Consider for instance a nonplanar double box in which the vertices that attach to external momenta k 1 and k 2 are both holomorphic. Then k 1 and k 2 are in fact orthogonal, The maximal cut of the all-massive nonplanar double box has two massless three-point vertices in a 2-mass-easy sub-box. Applying the above prescription suggests that we should expect four solutions, but in the previous section we found six. The explanation is that those solutions actually correspond to the four distinct distributions of chiralities, but two of them are in fact two-fold degenerate as is visible in fig. 4.
Once the pair of momenta in the sub-box contains at least one massless momentum, the counting works and the map is one-to-one as expected. In order to ensure healthy kinematics when more legs are massless, we add a supplementary condition: a valid diagram with nonchiral vertices is compatible with all massless limits for general momenta. The virtue of this rule is illustrated in fig. 6 for two graphs that do not appear in our hepta-cut solutions. To better understand this, suppose that the remaining leg in the sub-box in the left diagram is also massless. Then either k 1 · k 4 = 0 or k 3 · k 4 = 0 depending on the chirality of the fourth vertex. In practice, it is probably easiest to obtain the full variety of diagrams for a given configuration of massive and massless legs from the purely massless case [50] by introducing nonchiral vertices and eventually removing redundant graphs.

Residues of the Loop Integrand and Topological Structure
The internal momenta have eight degrees of freedom in strictly four dimensions and therefore we are left with one free parameter, after imposing the hepta-cut constraints. This implies that after promoting the original real slice contours of integration to seven-tori encircling the simultaneous zeros of the denominators in the crossed box integral, the integrand is localized onto the Riemann sphere associated with one of the hepta-cut equations found in the previous section. The maximal cut involves an inverse Jacobian that has multiple poles in the remaining variable z. These poles are known as composite leading singularities. For integrals with tensor numerators, additional poles in the integrand must be taken into account. Once the pole structure is properly understood, we can define the hepta-cut integral by a weighted residue expansion over a minimal set of all poles. We return to this part of the computation shortly.
In order to actually perform the hepta-cut contour integrals, it is necessary to break the automatic satisfaction of the on-shell constraints ℓ 2 1 = ℓ 2 2 = 0 which instead should be imposed by localization. To that end it is advantageous to introduce new parameters ζ 1 , ζ 2 ∈ C and two null-vectors η 1 , η 2 with the obvious properties, and then shift the loop momentum parametrization such that In general, this choice is tied to the specific on-shell solution in question. An appropriate choice in any case is to take the masslessly projected legs because it simplifies the calculation. For instance, let η 1 = k ♭ 2 and η 2 = k ♭ 3 . Then for all four-mass solutions S 1 , . . . , S 4 it is clear that and therefore we obtain the desired solution ζ 1 = ζ 2 = 0 since ξ 1 ξ ′ 1 = 0 = ξ 4 ξ ′ 4 for general momenta. The same arguments show that a valid choice for solutions S 5 and S 6 is We also need to include the Jacobian J L × J R due to coordinate transformation from loop four-momenta to parameter space. This can easily be done by Wick rotation to Euclidean spacetime where the volume of the 4-parallelotope spanned by a set of four vectors q µ i can be computed from the corresponding Gram determinant up to an overall sign, which we eventually determine numerically along with potential factors of i introduced by the analytic continuation, The set of variables parametrizing the two loop momenta ℓ 1 and ℓ 2 after removal of redundant degrees of freedom also depend on the branch in question. In S 1 and S 3 we keep the variables α = (ζ 1 , ξ ′ 1 , ξ 2 , ξ ′ 2 ) and β = (ζ 2 , ξ 3 , ξ ′ 3 , ξ 4 ) and thus fix ξ 1 = ξ ′ 4 = 1. The corresponding Jacobians to appear in the numerator are then Notice that these forms are not constant because the parametrization is not linear in the parameters as opposed to previous work [50] in the purely massless case. Similarly for S 2 and We finally compute the Jacobians appropriate to solutions S 5 and S 6 , (4.8)

Hextuply Pinched Genus-3 Curve
The preceding discussions now lead us to the definition of localization of the two-loop nonplanar double box integral onto the Riemann sphere associated with the ith on-shell branch. Changing integrations variables from loop momenta (ℓ µ 1 , ℓ µ 2 ) to parameters (ξ i , ξ ′ i ), replacing real slice integration contours by a multidimensional toris encircling the joint solution of the hepta-cut constraints and subsequently performing seven contour integrals using eq. (2.9) give rise to a total Jacobian J i , The generic form of this Jacobian is a product of n i simple-pole factors associated with the pinching points or intersections with neighbouring Riemann spheres, where h(z) is a regular function of z. For the two-loop crossed box integrals with up to four massive legs and no doubled propagators, the Jacobian will at most define a quartic polynomial because n i ≤ 4 for all i as we shall see below.
It is straightforward to obtain the Jacobians explicitly after the inverse propagators have been expanded in parameter space. For brevity we merely state the results here. We refer the reader to e.g. refs. [46,50,51] for related examples. In advance of calculations below, we identify a frequently occuring kinematical constant along the lines of ref. [49], (4.11) The multivariate residues evaluated at the simultaneous zeros of the denominators are, (4.13) (4.14) where the pole locations in the Jacobians are directly exposed, This leads to the definition of the octa-cut of a general tensor integral with numerator insertion Φ(ℓ 1 (z), ℓ 2 (z)), Here, Γ i is a weighted linear combination of small circles around the poles in the remaining variable z chosen so that the integrals extract the residues of the loop integrand. Denoting the weight of the residue evaluated at z = ξ for the ith branch by ω(i, ξ), we have Not all of these residues are independent though, as can be explained from the global structure of the unitarity cut [47,64]. Indeed, consider an arbitrary integrand test function of the two loop momenta, Φ(ℓ 1 (z), ℓ 2 (z)), and assume regularity on the Jacobian poles. It is then very easy to prove the residue relation, where the left and right hand sides of the equation are understood to be evaluated in local coordinates on solutions S i and S j respectively. Other choices are equally valid, e.g symmetric in i and j. For the purpose of completeness, let us state all such identities: This pattern of intersections confirms the global topological structure of the hepta-cut in fig. 8. This picture follows by pinching the tubes of the genus-3 surface six times along a horizontal and a vertical line passing through the center of the object, see fig. 7. The number of independent residues is reduced to 16. The residues at infinity and in numerator insertions are not shared. However, within each Riemann sphere the one-dimensional global residue theorem ensures that the sum of all residues vanish. 1 ∩ 4 and at least one massless particle (right). All degenerate configurations considered in this paper fall within these two topological pictures. The straight lines are drawn for simplicity and should be interpreted as Riemann spheres.
The upshot is that we only need to encircle 10 global poles to produce a complete basis of homology for S 1 ∪ · · · ∪ S 6 . In particular, the contributions from branches S 3 , . . . , S 6 are redundant, because each of them only has a single pole (at z = ∞) besides those located at intersections with S 1 and S 2 . To minimize an overcomplete basis we may simply set extraneous residue weights to zero. It is natural to choose an ordered set of winding numbers, call it Ω, so that we encircle all Jacobian poles along with poles in numerator insertions where both loop momenta become infinite simultaneously. Following the notation of ref. [47], the weights of the ten global poles can be written (4.20) By convention, a residue with weight ω i∩j is evaluated on the i'th branch. Later it may be convenient to instead encircle infinity poles, since scalar integrals have simpler residues there.

Octuply Pinched Genus-3 Curve
Let us now relax the condition m 1 m 2 = 0 and analyze the analytic structure of the loop integrand under those circumstances. The eight hepta-cut solutions for this class of kinematics were determined in section 3.3. The topological picture is that of an octuply pinched genus-3 surface, where tubes have been contract along one vertical and two horizontal lines through the center as shown in fig. 7. We can reproduce this situation locally from coincidence of residues. Localizing the scalar master integral onto each of the eight Riemann spheres paramatrized by the hepta-cut solutions S 1 through S 8 yields four pairs of one-dimensional contour integrals, , (4.23) The residues at the poles in the displayed loop integrands again satisfy a rich set of linear relations across the irreducible components of the genus-3 curve and thus reflect the global structure of the unitarity cut. In fact, where χ ≡ s 14 /s 12 . We point out that the hepta-cut contributions in eqs. (4.21)-(4.24) in ref. [50] after appropriate reparametrization, . (4.28)

Master Integral Projectors
The hepta-cut localizes the nonplanar double box integral onto a variety of linked Riemann spheres associated with the joint solutions of the on-shell equations. What remains is a one-dimensional complex contour integral whose integrand has poles. We can now choose contours that extract residues of the integrand and effectively obtain an octa-cut such that the integral is completely localized to a point in C 8 . On top of that, consistency of the unitarity method imposes nontrivial constraints on these contours, however. The reason is that converting a real slice integral into a multidimensional contour integral in general does not respect various relations among integrals. Each identity leads to a constraint which at two loops always can be rearranged and phrased as vanishing of a certain function upon integration over R D × R D . Schematically, It is easy to understand the nature of these relations if we imagine that we compute an amplitude diagram by diagram. All contractions between loop momenta and external vectors are expressible in terms of the eight fundamental scalar products ℓ i · e j for i = 1, 2 and j = 1, 2, 3, 4 where e = (k 1 , k 2 , k 4 , ω). Here, ω is a spurious direction that is perpendicular to the subspace spanned by the four external momenta. Odd powers of ℓ 1 ·ω and ℓ 2 ·ω vanish upon integration, whereas even powers are reducible in four dimensions. It readily follows that ℓ 1 · k 1 , ℓ 1 · k 2 and ℓ 2 · k 4 can be written in terms of inverse propagators and external invariants. Moreover, ℓ 2 · k 2 depends linearly on ℓ 1 · k 4 and ℓ 2 · k 1 . The latter two may be selected as irreducible scalar products. Accordingly, the general numerator polynomial for the problem at hand can be parametrized as follows, The integrand reduction can be obtained by imposing renormalizability conditions that constrain the exponents of the ISPs and performing the multivariate polynomial division of N modulo a Gröbner basis constructed from the seven inverse propagators. The latter part including identification of ISPs is carried out automatically by the Mathematica package BasisDet [57]. The integrand contains 19 parity-odd and 19 parity-even elements as previously reported [54]. The analysis above suggests that the amplitude contribution in question contains 19 genuine integrals of the form Many integrals are expressible as linear combinations of integrals with lower-rank tensors and fewer than seven propagators. This reduction is achieved due to IBP relations that follow from inserting a total derivative into the loop integrand, where · · · means fewer-propagator topologies that have vanishing hepta-cuts, and the consistency constraint thus reads The IBP relations can be generated by various public computer codes. For this project we used the Mathematica package FIRE [70]. We list below a few examples. There are two relations which hold for arbitrary values of the external masses, The parity-odd terms in the integrand basis of course vanish identically after the loop momentum integration has been performed, but they are nonetheless very important for integrand-level reduction and unitarity purposes. It is sufficient to require that the full variety of integrals with Levi-Civita insertions, after invoking momentum conservation, ε(ℓ 1 , k 2 , k 3 , k 4 ) , ε(ℓ 2 , k 2 , k 3 , k 4 ) , ε(ℓ 1 , ℓ 2 , k 1 , k 2 ) , ε(ℓ 1 , ℓ 2 , k 1 , k 3 ) , ε(ℓ 1 , ℓ 2 , k 2 , k 3 ) , (5.11) continue to integrate to zero on general contours in C 4 × C 4 . Our goal of the rest of this paper is to massage the amplitude master equation (2.1) into a form that allows us to project the master integral coefficients.

One-Mass Projectors
The simplest configuration is the one-mass diagram with, say, m 1 = 0. As in the purely massless case [50,54], there are two master integrals so the amplitude contribution can be expressed as where integrals with less than seven propagators are truncated. Since momentum k 4 is massless, we have ℓ 1 · k 4 = (ℓ 1 + k 4 ) 2 /2. We decide to encircle the following set of global poles, and let Ω denote the corresponding weights, . . . , ω 10 ) . (5.14) The hepta-cut two-loop crossed box integrals reduce to  We exploit the freedom to choose contours after imposing the reduction conditions and define two master integral projectors (also called master contours) which extract of either of the master integral coefficients, Here, M 1 and M 2 are just particular lists of the winding numbers of the corresponding global poles with the property that they only receive contribution from one master integral, which is also normalized to unity. The eight contour constraints together with either of the projectors are in practice arranged as 10 × 10 matrices. The rank is 10 and the solutions for the weights are uniquely determined. We find that the projectors are characterized by the 10-tuples The master integral coefficients can be written compactly in terms of tree-level data as where the rescaled Jacobian for this configuration is defined by J(z) ≡ ± m 2 1 − s 12 χs 12 1 z(z + 1)(z + (s 12 (1 + χ) − m 2 1 )/(χs 12 ) .

(5.27)
The computation of the remaining one-mass configurations is essentially equivalent to the one described here and the projectors are similar. The lack of symmetry in the nonplanar double box suggests that we also derive projectors for the one-mass diagram with m 4 = 0 (or m 3 = 0). The hepta-cuts in the limit m 1 , m 2 , m 3 → 0 follow from eqs. , (5.29) whereas the singular point locus and parity-odd vanishing constraints carry over direclty from the calculation above. The residues in the masters are  The expressions for the one-mass projectors derived here are consistent with the purely massless calculation reported in ref. [50].

Two-Mass Projectors
As previously explained, there are four kinematically inequivalent distributions of massless and massive external legs in the two-mass four-point nonplanar double box. Indeed, we distinguish between the two-mass short-side diagrams with either both massive legs situated in the planar or nonplanar end and the long-side and diagonal diagrams. From the point of view of the global structure of the hepta-cut, the latter three are similar and can be treated within the regime of the octuply pinched genus-3 curve whereas the first diagram is a variant of the three-and four-mass case.
The long-side two-mass diagram can be studied by taking over the singular point locus (5.14), basis integral decomposition as well as the parity-odd contour constraints (5.17) it turns out. We assume that m 1 m 4 = 0 and m 2 = m 3 = 0. Under these circumstances, the relevant hepta-cuts are , (5.38) where the pole location λ is defined by The residues associated with the two master integrals are as follows, where N 1 and N 2 are given by (q 1 , q 2 , −q 1 , −q 2 , q 3 , q 1 , q 2 , −q 1 , −q 2 , q 3 ) , (5.46) for constants q 1 , q 2 and q 3 where We have also derived projectors for the two-mass diagonal configuration. The computation essentially resembles that of the two-mass longside diagram, meaning that the same singular point locus, integral basis and contour constraints can be used. For m 2 = m 4 = 0 and m 1 m 3 = 0, the hepta-cuts evaluated at the branches S 1 , . . . , S 4 are , (5.54) where the third pole λ is now defined by It can be shown that the master contours for this kind of integrals are characterized by the following numbers, Our next example is the short-side two-mass diagram with m 1 m 2 = 0 and m 3 = m 4 = 0 which is a smooth limit of the three-and four-mass case. Accordingly, we now have the master equation and we must choose a new singular point locus for this computation. It is natural to encircle the following ten global poles, We remind that (G 1 , . . . , G 5 ) and (G 5 , . . . , G 10 ), are located at the following values of z on the Riemann sphere parametrizated by S 1 and S 2 respectively, Evidently, on-shell branches three through six are eliminated and we retain only the following two hepta-cut integrals, , (5.62) where the overall constant is given by  Resolving the contour constraints leaves three contour weights undetermined, exactly balancing the number of master integrals for this problem. The master contours which pick up contribution from one basis integral and annihilate the other two are, After changing the remaining one-dimensional contour into a linear combination of small circles around the global poles, the three master integrals reduce to the following residues, Res {G i } X * * 2,1,1 [1] = N 1 (r 1 , −r 2 , −r 1 , r 2 , 0, r 1 , −r 2 , −r 1 , r 2 , 0) , (5.69) where N 1 ≡ − m 2 1 m 2 2 + χγ 12 s 12 χs 12 (γ 12 + m 2 2 )(γ 2 12 + m 2 1 m 2 2 + χγ 12 s 12 ) , (5.72) , (5.73)  where q 1 , . . . , q 5 are defined by q 1 ≡ γ 2 12 − m 2 1 m 2 2 − γ 12 χs 12 , q 2 ≡ − 2m 2 1 m 2 2 − γ 12 χs 12 , q 3 ≡ γ 2 12 − m 2 1 m 2 2 + γ 12 χs 12 , q 4 ≡ q 1 − q 2 + q 3 q 5 ≡ q 1 + q 3 . (5.79) We point out that although the projectors are functions of irrational quantities such as γ 12 which has a square root, the final integral coefficients obtained in this way are rational in external invariants.

Three-Mass Projectors
The three-mass case with m 1 = 0 or m 2 = 0 and the remaining three external masses nonzero is similar to the two-mass long-side calculation previously presented. Here we will instead focus on the three-mass diagram with m 4 = 0 which has the two-mass short-side diagram with m 1 m 2 = 0 as a smooth limit. This means that we can continue to use the integral basis (5.59), the singular point locus (5.60) and the parity-odd constraints (5.65). Moreover, the contour constraints from IBP relations are identitical to those in eqs. (5.66)-(5.67).
For this problem, we will use the hepta-cuts for S 1 and S 2 from eq. (4.14) with µ = 0 and without the overall factor of γ ⋆ , It is possible to obtain rather clean forms of the residues computed by the hepta-cut master integrals at the singular point locus, if we prefer quantities constructed from flat momenta, ξ 1 ,ξ 2 and τ , instead of the usual Mandelstam variables and external masses. Indeed, the residues can be expressed as , , (5.83) The master contours that respect all integral reduction identities and extract either of the three master integrals are In these equations, the overall constants N 1 , N 2 and N 3 are , (5.89) along with the residue weights q 1 , . . . , q 5 , We note that the three-mass projectors written here reduce to the two-mass shortside formula (5.68) in the limit m 3 → 0.

Four-Mass Projectors
We finally examine the principal kinematical configuration with four distinct external masses m 2 i = k 2 i = 0. The intermediate calculations are more complicated because neither γ 12 nor γ 14 can be simplified to rational expressions. In order to simplify the computation to the maximum extent possible prior to solving for the projectors, we will encircle a slightly different set of global poles compared to the previous examples. More specifically, we exploit that integrand of the scalar master integral, evaluated on branches three through six, have vanishing residues at infinity.
It is convenient to arrange the poles so that G i and G i+5 are still parity conjugates of each other. The set of global poles is The residues at these poles can be streamlined by rescaling all hepta-cut Jacobians J 1 , . . . , J 6 by a common factor, and this is what we will do implicitly below. This constant will eventually drop out when we compute integral coefficients.
Without repeating the exercise, we know that the contour constraints from parity-odd numerator insertions are Prior to presenting the residues computed by the three master integrals, it proves advantageous to define the following four constants constructed out of various previously defined quantities, In terms of the r i 's, the residues of the master integrals can be brought to a particularly simple form, At first sight these constraints differ from those found in the two-and three-mass calculations, see e.g. eqs. (5.66)-(5.67). However, the two pairs of equations enforce the same constraints, as can be argued easily. We may express the constraints without imposing the global residue theorem. For the four-mass case we then have and Then it is immediately clear that we obtain the same answer from eqs. (5.107)-(5.108) when we truncate to linearly independent sets of residues by the global residue theorem.

Reduction of Integrals with Doubled Propagators
Feynman integrals with doubled and in general higher powers of propagators frequently appear in loop amplitude computations, for instance in IBP identities, Schwinger parametrizations or bubble insertions. It was recently explained that generalized unitarity cuts of such integrals are naturally treated as degenerate multivariate residues using computational algebraic geometry [52]. In that connection, several examples were given for one-and two-loop integrals with massless kinematics. This method extends seamlessly to multiloop integrals with external masses, as we will now demonstrate for the two-loop crossed box.

Unitarity Cut Algorithm: Bezoutian Matrix Method
We very briefly review the unitarity cut algorithm for integrals with higher powers of propagators. For more details and examples, please refer to ref. [51,52]. The main ingredient needed is computational algebraic geometry.
Recall that a residue is nondegenerate, if the Jacobian at the pole ξ is nonzero, i.e., In this case, the value of the residue is simply calculated by Cauchy's theorem in higher dimensions, i.e. eq. 2.9. However, the Jacobian clearly vanishes if there is one or more doubled propagators being cut and the residue is degenerate, so this approach does not apply. To solve the problem, we need techniques from algebraic geometry. There are two ways of evaluating such residues: 1. The transformation law (see for instance ref. [72]). This theorem can be used to convert a degenerate residue at the simultaneous zero of the inverse propagators to a factorizable residue. The explicit transformation matrix is found by the Gröbner bases method. The algorithm is described in refs. [51,52].
2. The Bezoutian matrix method. Here one determines the duality structure [71] of the multivariate residues, which in turn can be calculated easily. In general, the Bezoutian matrix method is considerably faster than the transformation law for complicated cuts with many independent external invariants.
Our Mathematica package MathematicaM2 2 is capable of computating multivariate residues using either of these techniques. In what follows, we outline the Bezoutian matrix approach and provide some basic examples.
Let I = f 1 , . . . , f n be an ideal in the ring R = C[z 1 , . . . , z n ]. Assume that I is a zero-dimensional ideal, i.e. the zero locus Z(I) = {ξ 1 , . . . , ξ k } consists of finite number of discrete points. For a zero-dimensional ideal I, the quotient ring R/I is a finite dimensional C-linear space.
Before we calculate individual residues, we first examine the structure of the sum of residues by Bezoutian Matrix. Then we eventually get individual residues from partition functions. For a polynomial h in R, we define the global residue as which is just the sum of all residues. By Stokes' theorem, the values of the residues only depend on h's equivalence class [h] in R/I. Furthermore, we can define an inner product , in R/I, g, h ≡ Res(g · h) . (6.3) Theorem 1 , is a nondegenerate inner product in R/I.
The proof of the theorem is given in ref. [72]. This theorem implies that, given a linear basis {p i } for R/I, we can find its dual basis {∆ i } in R/I, such that In practice, the basis and dual basis can be found by the Gröbner basis method and the Bezoutian matrix [74]. The procedure involves the following steps: 1. Calculate G, the Gröbner basis of I in the DegreeLexicographic order. Denote the leading terms for all polynomials in G as LT (G). Then all monomials in R which are lower than LT (G) constitute {p i }, which is the canonical linear basis for R/I.

2.
Introduce a set of auxiliary variables {y 1 , . . . , y n } and define the Bezoutian matrix B for I as, Calculate its determinant, det B. The dual basis explicitly characterizes the structure of global residues. Let the decomposition of the unit 1 over the dual basis be given as Then for an arbitrary numerator h, expand [h] over the canonical linear basis, 8) and the global residue is given as [74], This formula is the result of the definition of the dual basis, and provides a very efficient way of calculating the residues.
To get individual residues, we can use the formula (6.9) and the new ingredient partition functions of Z(I). This theorem can be proven by construction [75]. Then for each individual residue at ξ i [74], we have the result Explicitly, the partition functions e 1 , . . . , e k can be constructed by the method of Lagrange interpolation. The computation via Bezoutian matrix method is realized in our package, MathematicaM2. We demonstrate this computation by a simple example before we return to generalized unitarity cuts of integrals with doubled propagators.

Example: One-Mass Nonplanar Double Box
Let now us apply the Bezoutian matrix algorithm to generalized unitarity cuts of integrals with doubled propagators. We will be slightly more general than in eq. (3.2) and define the two-loop crossed box integral with arbitrary integer powers (σ 1 , . . . , σ 9 ) of propagators and irreducible numerators as where the seven propagators f 1 , . . . , f 7 can be found in eq. (3.1) with k 5 = k 6 = 0 and In order to unambigously define the degenerate multivariate residue associated with the maximal cut, the inverse propagators are grouped into seven factors g i ≡ f σ i i . As in [52] we will for technical simplicity only consider cuts in strictly four dimensions, postponing the analysis in D = 4 − 2ǫ dimensions to future work.
Our principal result is unique analytic contours for all basis integral coefficients in all inequivalent configurations of massive and massless external momenta in the four-point two-loop crossed box, valid to O(ǫ 0 ) in the dimensional regulator. The content of this paper is also relevant for higher-multiplicity scattering of massless particles. The maximal cut defines a nodal algebraic curve associated with a hextuply or octuply pinched genus-3 Riemann surface whose components are Riemann spheres. The first category includes contributions where both legs in the planar end of the diagram are massive, whereas the second covers the rest. The number of linearly independent residues is always ten as expected. We find that for the sixfold degenerate curve, the projectors for all three master integrals are unique once we impose five linearly independent Levi-Civita constraints and two linearly independent IBP conditions. The Levi-Civita constraints are resolved for weights that respect parity. In the four-mass case, unlike the situation for the four-mass double box [49], the IBP constraints are not satisfied automatically for the nonplanar double box.  Table 1. Classification of all kinematically disctinct diagrams from the viewpoint of the maximal cut. The columns list whether the external masses m 1 , m 2 , m 3 and m 4 are zero or not, the number |S| of hepta-cut solutions, the number of independent residues, the number of parity-odd and parity-even contour constraints and finally the set of master integrals. The notation for the master integrals refers to the powers of the two irreducible numerator insertions.
Overall, the results exhibit a very interesting and naively unexpected simplicity, which clearly deserves more attention. Indeed, the systematics of the contour constraints (e.g. eqs. (5.66)-(5.67)) are remarkable. Instead of being a set of disconnected calculations, the inequivalent kinematical configurations related through a rich underlying structure that seems to be governed by the global picture of the hepta-cut. As summarized in table 1, we find that the contour constraints are identical for all configurations within a particular class of hepta-cut solutions, for example all the way from the three-mass diagram with m 1 m 2 = 0 to the purely massless case. The IBP contours seem to be even more systematic. In all cases with six hepta-cut solutions we find two linearly independent IBP constraints. The chiral branching from 6 → 8 hepta-cut solutions triggers the emergence of an additional IBP constraint. Interestingly, as for the planar double box [49], the IBP constraints are inheritered through chiral branchings. To see this, let us instead consider hepta-cuts from the eightfold degenerate genus-3 curve, with the set of global poles, ( G 1∩6 , G 1∩8 , G 3∩5 , G 3∩7 , G 3,∞ R , G 2∩5 , G 2∩7 , G 4∩6 , G 4∩8 , G 4,∞ R ) . (7.1) As pointed out in appendix B, these poles exactly correspond to poles on the sixtuply pinched genus-3 curve. It is now an easy task to check that two of the three IBP constraints are inherited. The same observation applies to the four-mass computation.
In view of the complexity of the hepta-cut expressions and the typical amount of effort required to generate IBP relations for high-rank integrals with many external relations, it is also striking that the constraints coefficients are simply integers. This also applies to the planar double box with up to four external masses [46,48,49] and the planar triple box [50]. This is a clear hint of a general principle that may be explained by algebraic geometry.
The last part of this paper described a new algorithm based on the Beouzitan matrix and Gröbner bases to compute degenerate multivariate residues which typically appear in more complicated calculations. This algorithm was applied to the reduction of a massive nonplanar double integral with a doubled propagator onto master integrals with single propagators. We expect the Bezoutian method to become increasingly valuable for problems involving two-loop topologies with fewer propagators and at three loops and beyond (see e.g. ref. [51]).
We end this paper by suggesting interesting projects for future research. It is desirable to understand the nature of the contour constraints in complete detail. In particular, is it possible to fully determine constraints arising from integration-by-parts identities directly from the underlying algebraic geometry? Recent progress for the planar double box shows that discrete symmetries to some extent determine these constraints [49]. Such symmetries seem to be less constraining at higher genera. We also find it urgent to extend maximal unitarity to D dimensions to recover terms missed in strictly 4D. Another very important next step is to extend the method to basis integrals with five external legs, for example the pentabox and turtle-box and the related nonplanar diagrams. First of all from a phenomenological point of view, but we also hope that a generalization beyond four external particles will offer insight in uniqueness of projectors [47]. We expect that the onemass hepta-cuts presented in this paper to be valuable in that direction, because octa-cuts may be evaluated as hepta-cuts followed by a particular choice of contour that puts the last propagator on-shell. Theoretically speaking, the master integral coefficients for all-massive six-point planar and nonplanar double boxes are exciting to compute because in those cases the on-shell parametrization is irrational and hence the maximally cut integrals suffer from genuine branch cuts [47]. Ultimatively, it would be intriguing to implement the formalism numerically. We are looking forward to address some of these questions soon.
in draft stage and helpful comments. Both authors are grateful to Institut de Physique Théorique, CEA-Saclay, and especially David Kosower for hospitality during stages of this project. The work of YZ is supported by Danish Council for Independent Research (FNU) grant 11-107241. and similarly for the residues at z = 0 and z = ∞, (G 11 , G 12 , G 5 , G 10 , G 13 , G 14 , G 15 , G 16 )ξ 2 →0 − −− → ( G 11 , G 12 , G 13 , G 14 , G 15 , G 16 , G 17 , G 18 ) . (B.6) The remaining two Jacobian global poles G 5 and G 10 are located at the nodal points S 1 ∩ S 3 and S 2 ∩ S 4 respectively, and are thus generated by chiral branching from 6 → 8 hepta-cut solutions. By the one-dimensional Global Residue Theorem, the number of independent residues is clearly invariant.