Maximal Unitarity for the Four-Mass Double Box

We extend the maximal-unitarity formalism at two loops to double-box integrals with four massive external legs. These are relevant for higher-point processes, as well as for heavy vector rescattering, VV ->VV. In this formalism, the two-loop amplitude is expanded over a basis of integrals. We obtain formulas for the coefficients of the double-box integrals, expressing them as products of tree-level amplitudes integrated over specific complex multidimensional contours. The contours are subject to the consistency condition that integrals over them annihilate any integrand whose integral over real Minkowski space vanishes. These include integrals over parity-odd integrands and total derivatives arising from integration-by-parts (IBP) identities. We find that, unlike the zero- through three-mass cases, the IBP identities impose no constraints on the contours in the four-mass case. We also discuss the algebraic varieties connected with various double-box integrals, and show how discrete symmetries of these varieties largely determine the constraints.


Introduction
Last year's discovery [1,2] by the ATLAS and CMS collaborations of a Higgs-like boson completes the particle content of the Standard Model. Coupled with the absence to date of direct signals of physics beyond the Standard Model, the discovery points towards an important role for precision measurements in determining the scale of new physics beyond the Standard Model.
Theoretical calculations at the LHC, whether for signals or backgrounds, begin with the tree-level amplitudes required for leading-order (LO) calculations in perturbative quantum chromodynamics (QCD). Because the strong coupling α s is relatively large and runs quickly, LO predictions suffer from strong dependence on the unphysical renormalization and factorization scales and are thus not quantitatively reliable. Next-to-leading order (NLO) is the lowest order in perturbation theory which offers quantitatively reliable predictions. These calculations require one-loop amplitudes in addition to tree-level amplitudes with higher multiplicity. Recent years have seen major advances in NLO calculations, especially for processes with several jets in the final state [3][4][5][6][7][8][9]. While the uncertainty left by scale variation cannot be quantified in the same fashion as statistical uncertainties, experience shows that it is of O(10-15%).
As combined experimental uncertainties in future measurements push below this level, a comparison with theoretical calculations will require pushing on to next-to-next-to-leading order (NNLO) accuracy. Such studies will require computation of two-loop amplitudes. These computations form the next frontier of precision QCD calculations. The only existing fully-exclusive NNLO jet calculations to date are for three-jet production in electron-positron annihilation [10]. These calculations have been used to determine α s to 1% accuracy from jet data at LEP [11]. This extraction is competitive with other determinations. Beyond their use in seeking deviations in precision experimental data from Standard-Model predictions, NNLO calculations will also be useful at the LHC for improving our understanding of scale stability in multi-scale processes such as W +multi-jet production, as well as for providing honest theoretical uncertainty estimates for NLO calculations.
The unitarity method [12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29] has made many previously-inaccessible one-loop calculations feasible. Of particular note are processes with many partons in the final state. The most recent development, applying generalized unitarity, allows the method to be applied either analytically or purely numerically [30][31][32][33][34][35][36][37][38][39][40]. The numerical formalisms underlie recent software libraries and programs that have been applied to LHC phenomenology. In this approach, the one-loop amplitude in QCD is written as a sum over a set of basis integrals, with coefficients that are rational in external spinors, (1.1) The integral basis for amplitudes with massless internal lines contains box, triangle, and bubble integrals in addition to purely rational terms (dropping all terms of O(ǫ) in the dimensional regulator). The coefficients are calculated from products of tree amplitudes, typically by performing contour integrals via discrete Fourier projection. In the Ossola-Papadopoulos-Pittau (OPP) approach [21], this decomposition is carried out at the integrand level rather than at the level of integrated expressions.
Higher-loop amplitudes can also be written in the form given in eq. (1.1). As at one loop, one can carry out such a decomposition at the level of the integrand. This generalization of the OPP approach has been pursued by Mastrolia and Ossola [41] and collaborators, and also by Badger, Frellesvig, and Zhang [42]. The reader should consult refs. [43][44][45][46][47][48][49] for further developments within this approach. Arkani-Hamed and collaborators have developed an integrand-level approach [50][51][52][53][54][55] specialized to planar contributions to the N = 4 supersymmetric theory, but to all loop orders.
Within the unitarity method applied at the level of integrated expressions, one can distinguish two basic approaches. In a 'minimal' application of generalized unitarity, used in a number of prior applications and currently pursued by Feng and Huang [56], one cuts just enough propagators to break apart a higher-loop amplitude into a product of disconnected tree amplitudes. Each cut is then a product of tree amplitudes, but because not all possible propagators are cut, each generalized cut will correspond to several integrals, and algebra will be required to isolate specific integrals and their coefficients. This approach does not require a predetermined general basis of integrals; it can be determined in the course of a specific calculation. A number of calculations have been done this way, primarily in the N = 4 supersymmetric gauge theory [57][58][59][60][61][62][63][64], but including several four-point calculations in QCD and supersymmetric theories with less-thanmaximal supersymmetry [65][66][67][68][69][70][71]. Furthermore, a number of recent multi-loop calculations in maximally supersymmetric gauge and gravity theories have used maximal cuts [72][73][74][75][76][77][78], without complete localization of integrands.
We will use a more intensive form or 'maximal' form of generalized unitarity. In this approach, one cuts as many propagators as possible, and further seeks to fully localize integrands onto global poles to the extent possible. In principle, this allows one to isolate individual integrals on the right-hand side of the higher-loop analog of eq. (1.1). In previous papers [79,80], we showed how to extract the coefficients of double-box master integrals using multidimensional contours around global poles. In this paper, we recast this operation as applying generalized discontinuity operators (GDOs). Each GDO corresponds to integrating the integrand of an amplitude or an integral along a specified linear combination of multidimensional contours around global poles. The GDOs generalize the operation of cutting via the Cutkosky-rule replacement of propagators by on-shell delta functions.
Some of the contour integrations in a GDO put internal lines on shell, equivalent to cutting propagators [81]. This integration will typically yield a Jacobian giving rise to poles in the remaining degrees of freedom. In the case of the double box, the Jacobian allows one to fully localize the remaining degrees of freedom through additional multidimensional contour integrals. The integrand is then fully localized at one of a set of global poles. We call these additional degrees of freedom 'localization variables'. We include these additional dimensions of contours in the definition of the GDO. This maximal-unitarity approach may be viewed as a generalization to two loops of the work of Britto, Cachazo, and Feng [16], and of Forde [23].
The GDOs are constructed so that each one selects a specific master integral, Applying it to eq. (1.1) then gives us an expression for the corresponding coefficent, The right-hand side will have the form of explicit contour integrals over localization variables of a product of tree amplitudes; schematically, The weights with which the contours Γ j surround the different global poles are determined by a set of consistency equations. These equations require that integrals vanishing over the Minkowski slice of complexified loop-momentum space are also annihiliated by the GDOs, vanishing on the particular combinations of contours in each and every GDO.
In this paper, we continue the maximal-unitarity approach of refs. [79,80], relying on the global-pole analysis [82] of Caron-Huot and one of the present authors. At higher loops, the coefficients of the basis integrals are no longer rational functions of the external spinors alone, but will in general depend explicitly on the dimensional regulator ǫ. We consider GDOs operating only on the four-dimensional components of the loop momenta, and accordingly extract only the leading terms, ǫ-independent terms. GDOs operating on the full D-dimensional loop momenta would be required to extract the remaining terms, and could presumably be used to obtain the rational terms in eq. (1.1) as well.
At two loops and beyond, the number of master integrals for a given topology will depend on the number and arrangement of external masses [83]. (See also recent work on a different organization of higher-loop integrals [84][85][86].) In previous papers, we have considered double boxes with no external masses [79] or with one, two, or three external masses [80]. In this article, we extend the GDO construction to planar double boxes with four external masses. We consider both the general case with unequal masses, and one special case with pairs of equal masses. In the general case, there are four master integrals; in the special case with an extra reflection symmetry, three. Søgaard [87] has constructed GDOs for the non-planar massless double box.
As in previous work, we ensure the consistency of the GDOs by requiring that they yield a vanishing result when applied to vanishing integrals. For the four-mass double box, it turns out that non-trivial constraints arise only from parity-odd integrands; integration-by-parts (IBP) identities [88][89][90][91][92][93][94][95] give no additional constraints. The symmetry requirement for the special equal-mass case must also be imposed explicitly, and unlike fewer-mass cases, does not emerge automatically from IBP equations. We consider only two-loop master integrals with massless internal lines. We will not consider the generalization to massive internal lines; but so long as there are sufficient massless internal lines to have at least one chiral vertex, the integrand should still have global poles, and we should expect the approach described here to generalize smoothly. This paper is organized as follows. In sec. 2, we present the parametrization of loop momenta we use for derivations. In sec. 3, we discuss the maximal-cut equations for the four-mass double box along with the global poles, and derive the GDOs for the four master integrals. In sec. 4, we discuss constraint equations and their symmetries for all double boxes from an algebraic-geometry point of view. We make some concluding remarks in sec. 5.

Loop-Momentum Parametrization
We take over the same loop-momentum parametrization used in ref. [80]. This parametrization makes use of spinors defined for massive external legs. Such spinors correspond to massless fourdimensional momenta, which we obtain using 'mutually-projected' kinematics. This construction was previously used in the work of OPP [21] and Forde [23] to extract triangle and bubble coefficients at one loop.
For a given pair of external four-momenta (k i , k j ), we require the mutually-projected momenta to satisfy (2.1) Figure 1: The double-box integral.
By construction, k ♭ i and k ♭ j are massless momenta. Next, define We note that so that ρ ij = k 2 i /γ ij . After using eq. (2.1), we obtain a quadratic equation for γ ij ; its two solutions are If either momentum in the pair (k i , k j ) is massless, only one solution survives. Eq. (2.3) then gives us γ ij = 2k i · k j . Inverting eq. (2.1), we obtain the massless momenta swap i ↔ j to obtain k ♭,µ j . In this paper we work with two mutually-projected pairs: (k 1 , k 2 ) and (k 3 , k 4 ). This choice defines a set of 'projected' massless momenta k ♭,µ i , i = 1, . . . , 4, in terms of the external momenta, k i , and the sign choices in γ ± 12 and γ ± 34 . With the projected momenta, we adopt the following parametrization for the on-shell doublebox loop momenta as depicted in fig. 1, This form ensures that ℓ 2 1 = ℓ 2 2 = 0. We write the various loop spinors in terms of the spinors corresponding to (k ♭ 1 , k ♭ 2 ) for ℓ 1 , and the spinors corresponding to (k ♭ 3 , k ♭ 4 ) for ℓ 2 : where the external spinors are defined via k ♭,µ i = i ♭ | γ µ |i ♭ ] /2. Without loss of generality, we can set two of the complex parameters to unity, ξ 1 = ξ 4 = 1, as we will do throughout the paper.
Moreover, similarly to ref. [80], we define the following quantities: where m i are the masses of the external momenta, m 2 i = k 2 i . In addition, we will make use of the following quantities not needed in ref. [80]: . (2.13)

Maximal Cuts of Double-Box Integrals
Our aim is to determine the coefficients of the double-box master integrals that appear in the basis expansion (1.1) of a two-loop quantity that may either be an amplitude, form factor, or correlator. Without loss of generality, we refer to the two-loop quantity as an amplitude. The double-box integral topology is illustrated in fig. 1, and defines the internal momenta p j . The integral is defined in dimensional regularization with D = 4 − 2ǫ as, where Φ denotes an arbitrary polynomial in the external and internal momenta. We refer to it as a numerator insertion. At one loop, all numerator insertions can be expressed as linear combinations of propagator denominators, external invariants, and parity-odd functions which vanish upon integration; but this is no longer true at two loops and beyond. At higher loops, some polynomials Φ are irreducible. Integrals with certain irreducible-numerator insertions can be related to others using IBP identities, but in general several will remain as master integrals. We seek formulas for the double-box coefficients to leading order in the dimensional regulator ǫ in terms of purely tree-level input. We begin by cutting all double-box propagators on both sides of eq. (1.1). This immediately eliminates all integrals with fewer than seven propagators, or with a different topology, as cutting an absent propagator yields zero.
Heuristically, we may imagine using the Cutkosky rules, and simply replacing the cut propagators by on-shell delta functions. On the left-hand side of the equation, we would then obtain, where A tree (v) denote the tree processes at each of the six vertices of the diagram in fig. 1, and p j denote the momenta flowing through each of the propagators. The cuts have also eliminated any potential infrared divergences, so we can take the four-dimensional limit for the integrand. On the right-hand side of eq. (1.1), we would obtain a sum over expressions of the form, If we interpret the expressions in eqs. (3.2) and (3.3) literally, however, we face a problem. The integrations in these equations receive contributions only from regions of integration space where the loop momenta solve the joint on-shell constraints, For generic external momenta, the solutions to these equations are complex. So long as the integrations are over real momenta (R 4 × R 4 ), we simply get zero. Equating the two expressions will yield 0 = 0, which is true but useless for extracting the coefficient in eq. (3.3). Instead of thinking of the loop integrals as integrals over real momenta, we can choose to think of them as integrals in complex momenta, ℓ i ∈ C 4 , taken along contours comprising the real slice, Im ℓ µ i = 0. Changing the contour then gives us an alternative way of imposing a delta-function constraint, one that is valid for complex as well as for real solutions.
The utility of reinterpreting delta functions as contour integrals was previously observed in the context of twistor-string amplitudes [96,97], and is also standard in more formal twistor-space expressions [98]. In one dimension, we seek to localize an integral, even if q 0 becomes complex. Cauchy's residue theorem gives us precisely such a localization if we replace δ(q − q 0 ) by 1 2πi 1 q−q 0 , and take the integral to be a contour integral along a small circle centered at q 0 in the complex q-plane. Analogously, a product of delta functions can be defined as a multidimensional contour integral, where the contour T ε (q 0 ) is now a torus encircling the simultaneous solution of denominator equations. For the simple form of the denominator here, the contour will be a product of n small circles (ε ≪ 1), T ε (q 0 ) = C ε (q 01 ) × · · · × C ε (q 0n ), each centered at q 0j . The simultaneous solution of the denominator equations is called a global pole. The question of what it means for a torus to encircle a global pole is much more subtle in higher dimensions than for a contour to encircle a point in one complex dimension; but the subtleties will play no role in the present article. There is one important respect in which the multidimensional contour integrals behave differently from integrals over delta functions, namely the transformation formula for changing variables. Given a holomorphic function f = (f 1 , . . . , f n ) : C n → C n with an isolated zero 1 at a ∈ C n , the residue at a is computed by performing the integral over a toroidal contour, whose general definition is T ε (a) = {z ∈ C n : |f i (z)| = ε i , i = 1, . . . , n}. This contour integral satisfies the transformation formula Unlike the conventional formula for a multidimensional real integral over delta functions, it does not involve taking the absolute value of the inverse Jacobian. This ensures that this factor is analytic in any remaining variables on which it depends, so that further contour integrations can be carried out. We use multidimensional contour integrals to define generalized discontinuity operators. The GDOs for the double box will be eightfold integrals taken over contours that are linear combinations of basis contours. Each basis contour encircles a single global pole, and we will refer to global poles and their encircling contours interchangeably. Applying a GDO means changing the contour of the integration from one over the real slice of C 4 × C 4 to one over the GDO's associated contour. We want seven of the eight contour integrations to correspond to the seven on-shell constraints p 2 j = 0; to do so, the contours must ultimately encircle solutions to these constraints. The integrands in eqs. (3.2) and (3.3) are left unchanged. Imposing the seven constraints leaves one complex degree of freedom. The heptacut constraints thus define a Riemann surface in C 4 × C 4 . As we will see below, this Riemann surface contains a number of poles. Their presence will allow us to freeze the remaining degree of freedom, by choosing an appropriate contour of integration for the corresponding localization variable. Before discussing the poles, however, we first review the structure of the Riemann surface. As discussed in ref. [82], the maximal-cut Riemann surface for the double-box integral is a pinched torus, with the number of pinches equal to twice the number of vertical double-box rungs that attach to an on-shell massless three-point vertex. An on-shell massless three-point vertex is either chiral or anti-chiral, enforcing a two-fold branching of the kinematical parametrization, and implying a pinching of the parameter space. A more careful analysis of the kinematical solutions shows that chiral vertices are (anti-)correlated across the vertical rungs of the doublebox integral. Hence, we classify the different types of pinches by their effect on the vertical rungs.

Kinematical Solutions, Jacobians and Global Poles
In previous work [80], we assigned double boxes to one of three classes (a), (b), or (c), according to whether an on-shell massless three-point vertex is connected to: (a) the middle rung, (b) the middle rung and one outer rung, (c) all three vertical rungs 2 . We treated the two latter classes in ref. [80]. Here, we consider class (a), corresponding to the four-mass double box, illustrated in fig. 1. In this class, the solutions to the heptacut equations (3.4) form a doublypinched torus, shown in fig. 2.
Each lobe of the doubly-pinched torus corresponds to one of two kinematical solutions, S 1 and S 2 . In terms of the loop momentum parametrization of eq. (2.7), both solutions have ξ ′ 1 =ξ ′ 1 , ξ ′ 4 =ξ ′ 4 , ξ 1 = ξ 4 = 1, and the remaining four variables (ξ 2 , ξ ′ 2 , ξ 3 , ξ ′ 3 ) take on the following values, where theξ ′ i are defined in eq. (2.9). The Jacobian that arises from changing the integration variables of eq. (3.1) to the ξ i , ξ ′ i and subsequently performing seven contour integrals, d 8 ℓ 7 j=1 1/p 2 j =⇒ dz J, takes the generic form where (z i,1 , z i,2 ) are the local coordinates of the intersection with the neighboring solution(s), More generally, the Jacobian evaluated on a Riemann sphere will always be a product of simplepole factors associated with the pinching points (also known as nodal points) on the sphere. As mentioned above, the heptacut of the double-box integral arises from performing seven of the eight contour integrals, and yields a Riemann surface given by the solution to the joint on-shell constraints (3.4). We are left with a single complex degree of freedom (or localization variable) z, and the freedom to choose a contour for its integration. In order to localize the integrand completely, we should have this last contour encircle a pole in z. As in classes (b) and (c) treated in ref. [80], such poles can arise from two sources: the Jacobian factor (3.9), or from the numerator insertions Φ in eq. (3.1), which introduce an additional dependence on z.
The Jacobian poles are the pinching points G 9,10 in fig. 2. Because these points are shared between different on-shell solutions, one must decide on a convention for the sphere on which the corresponding residue is to be evaluated. We adopt the convention of computing the residue on the sphere located on the anti-clockwise side of the Riemann surface. In fig. 2, for example, the residue at G 9 should be evaluated on S 1 ; and the residue at G 10 on S 2 . Furthermore, we choose the orientations on each Riemann sphere such that for any global pole G k ∈ S i ∩ S j , the residues evaluated on spheres S i and S j are equal in magnitude but opposite in sign. That is, for an arbitrary function f of the loop momenta one has, , (3.11) in agreement with the conventions of refs. [80,82]. Other choices of conventions are possible, but all will lead to the same final expressions for the two-loop integral coefficients.
The class (a) Jacobian takes the form , for solution S 2 . (3.13) We have rescaled the original Jacobian J (0) into J by analogy with eq. (4.6) of ref. [80]. Note that z + z − =ξ ′ 2 /ξ ′ 1 . In addition to the Jacobian poles, we may choose the contour for the remaining post-heptacut degree of freedom to encircle any of the points on the Riemann surface where a loop momentum becomes infinite. At such points, numerator insertions Φ(p i ) have a pole in z. We will refer to these points, shown as punctures on the spheres in fig. 2, as insertion poles. There are eight such poles, so that altogther we have ten global poles which the z contour integral may encircle, located at the following values of (ξ 2 ,ξ ′ 2 , ξ 3 ,ξ ′ 3 ): (3.14) In the above labeling, the poles (G 2j−1 , G 2j ), j = 1, . . . , 7 form parity-conjugate pairs. Because parity amounts to swapping chiralities • ←→ • , thereby rotating fig. 2 by an angle π, parityconjugate pairs always appear antipodally in the figure. We note that at the pinching points G 9,10 , the loop momentum flowing through the middle rung of the double box becomes soft, p 4 → 0. At the remaining global poles, either the left or right loop momentum goes to infinity in a particular direction. (See the appendix for a more detailed discussion.) Let us denote a contour consisting of a small circle around G j by C j . The set of circles around all of the ten poles in fig. 2 forms an overcomplete basis of contours for GDOs, and equivalently an overcomplete basis for homology. On each sphere we can use the fact that all residues sum to zero to eliminate any one contour C j in favor of the remaining ones. This is not sufficient, as we must impose additional consistency constraints on the linear combination of contours by which every GDO acts. We discuss these below. Retaining all contours instead of choosing a linearlyindependent subset does have the advantage of making manifest certain discrete symmetries, clarifying the structure of the additional consistency constraints. We examine this issue in more detail in sec. 4.
The truncation to a linearly independent homology basis can be achieved simply by setting the coefficients of certain contours to zero in every GDO. Not all truncations will lead to a valid basis, however; a basis must necessarily contain a contour encircling at least one of the pinching points G 9,10 . To understand why, consider the sum of all residues on S 1 plus the sum of all residues on S 2 . Both sums are zero by Cauchy's theorem, and this sum is therefore zero. On the other hand, the sum equals that over the insertion poles alone, 8 i=1 Res G i , because the contributions from the pinching points cancel owing to eq. (3.11). We thus conclude that the residues at the insertion poles always sum to zero, and the set of contours encircling insertion poles alone does not constitute a complete homology basis on S 1 ∪ S 2 .

Master Contours -General Four-Mass Kinematics
Generalized discontinuity operators for the planar double box are given as eightfold contour integrals, which factor into a sevenfold contour integral localizing the integrand onto the heptacut solution surface -the joint solution of the on-shell equations for all seven propagator momenta. The last contour integral is now a contour integral on that Riemann surface. The contour cannot be chosen arbitrarily, however. It is subject to the consistency requirement that it yield a vanishing integration for any function that integrates to zero on the original contour of integration for the Feynman integral, R D × R D . This ensures that two integrals which are equal, for example by virtue of non-trivial integral relations, have the same generalized cut, Examples of terms which integrate to zero on R D ×R D include parity-odd terms and total derivatives used in the integration-by-parts identities to reexpress a large set of formally-irreducible integrals in terms of linearly independent master integrals. We may write a general contour for a GDO as follows, where the ω i are complex coefficients, and where the sum is taken over a linearly independent homology basis (or over an overcomplete one). For double-box integrals belonging to class (a), it turns out that consistency with IBP relations imposes no constraints on the contour, and hence no constraints on the ω i . On the other hand, the vanishing integration of (parity-odd) Levi-Civita numerator insertions -such as ε(ℓ 1 , k 1 , k 2 , k 4 ) -results in the following constraints on the coefficients ω i : 2ω 1 − 2ω 2 − ω 9 + ω 10 = 0 2ω 3 − 2ω 4 − ω 9 + ω 10 = 0 2ω 5 − 2ω 6 − ω 9 + ω 10 = 0 2ω 7 − 2ω 8 − ω 9 + ω 10 = 0 .

(3.29)
One is not obliged to make use of the integral identity (3.24) and enforce the ensuing contour constraint (3.28); one could equally well expand the equal-mass amplitude in terms of the slightly overcomplete basis in eq. (3.18), with the associated master contours given in eq. (3.22). Indeed, since the energies of heavy particles follow a Breit-Wigner distribution, an amplitude involving four massive vector bosons (e.g., W Z → W Z) will typically be required only for unequal masses; only when taking the on-shell approximation would the equal-mass case arise.

Varieties Arising from Feynman Graphs
In this section we discuss the heptacut of the planar double-box integral, putting some of the observations in sec. 3 into the broader context of algebraic geometry. On-shell constraints are polynomial equations. Accordingly, their simultaneous solution defines an algebraic variety. Ref. [82] observed that the variety corresponding to setting all seven propagator momenta of the planar double box on-shell is a pinched torus, with the number of pinches equal to twice the number of double-box rungs that end on at least one three-point vertex. As mentioned in the previous section and in ref. [80], we denote integrals having one, The components of the pinched tori are Riemann spheres. These spheres are associated with distinct solutions to the joint on-shell constraints (3.4) and are characterized by the distribution of chiralities ( • or • ) at the vertices of the double-box graph. The fact that the number of pinches is always even is a reflection of the fact that the on-shell solutions always come in parity-conjugate pairs. At a pinching point, there is exactly one double-box rung whose momentum becomes collinear with the massless external momenta connected to the rung. For the original uncut double-box integral, such regions of the loop momentum integration typically produce infrared divergences, and the pinches can therefore roughly be thought of as remnants of the original IR divergences. In addition, the pinched tori contain a number of insertion points (for example, in fig. 3(a), the points G 1 , . . . , G 8 ) where one of the loop momenta becomes infinite. Because the order of the pole is related to the ultraviolet power counting of underlying integrals in the theory (taking into account fermi-bose cancellations), these insertion points can be associated, roughly speaking, with UV divergences.
The pattern of global poles in classes (a)-(c) can be understood as follows. Starting from fig. 3(a), we can imagine taking massless limits of external momenta, at each step having exactly one additional double-box rung end on a three-point vertex. Geometrically, each step adds a pair of pinches (nodal points)-the first step producing fig. 3(b), and the second producing fig. 3(c). Each pinch preserves the global poles already present, and adds a global pole at the location of the pinching point. Nonetheless, pinching leaves the number of independent global poles constant: while it creates a new global pole, it also creates a new Riemann sphere and hence adds a global residue constraint which allows one global pole to be eliminated. There are eight independent global poles in class (a), and the number remains eight in classes (b) and (c).
To be more specific, class (b) contains 12 global poles, as illustrated in fig. 3(b). The poles G 1 , . . . , G 10 are obtained by taking the limitξ ′ 3 → 0 (corresponding to either m 3 → 0 or m 4 → 0) of the class (a) poles in eq. (3.14). 4 As evaluating the limit of G 9 and G 10 is slightly subtle, we quote the result here: where ω 11,12,13,14 → 0 in class (a), and ω 13,14 → 0 in class (b). If we choose a homology basis consisting of parity-conjugate pairs of poles, these constraints are expressed by the simple geometric statement that a valid contour must be invariant under a rotation through π radians of figs. 3(a), 3(b), and 3(c), respectively. Consistency with IBP identities imposes less transparent constraints on maximal-cut contours. In class (b), there is a single IBP constraint which takes the form, whereas in class (c) there are two IBP constraints which take the form, ) The first class (c) constraint (4.6) is identical to the class (b) one in eq. (4.5). This suggests that these constraints arise during the pinchings that carry the doubly pinched torus depicted in fig. 3(a) into the quadruply pinched torus of fig. 3(b) and thence into the sextuply pinched torus of fig. 3(c). The transition from fig. 3(a) into fig. 3(b) involves two (parity-conjugate) pinches which one might at first expect to produce two constraints. However, as a valid contour must be parity-symmetric (4.4), we should really expect one independent constraint to arise from a double pinching. This constraint is accompanied by a second constraint arising from the double pinching that turns fig. 3(b) into fig. 3(c). This pattern offers hope that it may be possible to derive the IBP constraints (4.5)-(4.7) directly from the underlying algebraic geometry. Expressing the IBP constraints in an overcomplete basis of homology makes it clear that they cannot be determined from algebraic topology alone. For example, on the sphere S 4 in fig. 3(b), the poles G 3 , G 5 , G 7 may be freely relabeled among each other without changing the topology. In contrast, eq. (4.5) does not have this relabeling symmetry.

Discrete Symmetries of IBP Constraints
We observe that the class (b) IBP constraint (4.5) is symmetric under reflection of fig. 3(b) in the vertical axis passing through the poles G 1 and G 2 . More explicitly, eq. (4.5) is symmetric under the interchanges 5 The pattern of relative minuses in eq. (4.8) owes to the fact that, in our orientation conventions, the reflection flips the orientation of the pinching or 'IR' cycles, but preserves that of the insertion or 'UV' cycles.
Conversely, assuming the symmetry (4.8), one might ask to what extent it determines the IBP constraint. The most general IBP constraint invariant under eq. (4.8) takes the form a 1 ω 1 + a 2 ω 2 + a 3 (ω 3 + ω 8 ) + a 4 (ω 4 + ω 7 ) + a 5 (ω 5 + ω 6 ) + a 6 (ω 9 − ω 11 ) + a 7 (ω 10 − ω 12 ) = 0 . (4.9) For convenience, let us now choose a basis of homology, for example ω 1,2,5,6 = 0. In this basis, the IBP constraint (4.9) takes the form, where we furthermore imposed the Levi-Civita constraints (4.4). Remarkably, the only thing left unexplained by the flip symmetry (4.8) is the fact that r Out of the two IBP constraints in class (c), we observe that eq. (4.6) is inherited directly from eq. (4.5) whereas the difference between eq. (4.6) and eq. (4.7) is symmetric under reflection of fig. 3(c) in a line passing through the centers of the spheres S 5 and S 6 . More explicitly, the difference is symmetric under the interchanges This symmetry will be broken by any choice of homology basis, highlighting the virtue of expressing the IBP constraints (4.6)-(4.7) in an overcomplete basis.

Conclusions
In this paper we have extended the maximal-unitarity formalism at two loops to double-box integrals with four massive external legs. We have constructed generalized discontinuity operators which isolate each of the four master integrals, annihilating all others. Applying one of these GDOs to the amplitude yields a formula for the corresponding coefficient in eq. (1.1), as a contour integral over products of tree-level amplitudes. We can choose to think of each GDO as operating in two steps. In the first step, we perform seven of the eight contour integrals, thereby putting on shell all internal lines of the double-box integral. This restricts the integrand to a Riemann surface, which has the form of a multiplypinched torus. Each component is a Riemann sphere, with the number of spheres equal to twice the number of double-box rungs that end on a three-point vertex. This step is identical for all four GDOs.
In the second step, we perform the remaining contour integral over a contour on the Riemann surface. This fully localizes the integrand onto a combination of global poles. The integration contours are different for each GDO. They are subject to consistency constraints. These constraints fall into two classes for general double-box integrals: (a) parity symmetry, amounting to invariance of the contours under rotations through π radians of the pinched tori; and (b) consistency with IBP relations. Writing out the latter constraints in an overcomplete basis of homology exposes additional flip symmetries. These symmetries alone would allow us to determine the constraints up to a small number of constants. The IBP relations determine these constants.
For the four-mass double box (class (a)), the underlying Riemann surface consists of two spheres, and there is no contour constraint from IBP relations. For the three-mass and short-side two-mass double boxes (class (b)), considered previously in ref. [80], the Riemann surface consists of four linked spheres. It can be viewed as derived from the two-sphere Riemann surface via a double pinching. One IBP constraint arises here. This constraint is inherited by the last case, a six-sphere surface corresponding to massless, one-mass, diagonal and long-side two-mass double boxes (class (c)), considered previously in refs. [79,80]. The six spheres again can be viewed as derived from the four-sphere surface via a double pinching, and an additional IBP constraint emerges as well. Thus, the IBP contour constraints appear to arise during the chiral branchings of the on-shell solutions, suggesting a strong connection to the underlying algebraic geometry.
This sequence of IBP constraints suggests a more natural choice of master integrals than that of eq. (3.18), namely double-box integrals with numerator insertions, 1, ℓ 1 · k 4 , ℓ 1 · k 4 − ℓ 2 · k 1 , (ℓ 1 · k 4 )(ℓ 2 · k 1 ) − 3 4 s 12 (ℓ 1 · k 4 ) − 1 8 s 12 s 14 . (5.1) These are linearly independent in class (a); in class (b), the last numerator insertion yields zero after integration, by virtue of an integration-by-parts relation; similarly, in class (c), the last two insertions produce vanishing integrals. A complete calculation of four-point amplitudes will also require the O(ǫ) terms in integral coefficients, and also GDOs for integrals with fewer than seven propagators. For processes with additional external legs, higher-point integrals will be needed as well. The simplest extension would probably be to 'turtle-box' integrals (P * 2,2 in the notation of ref. [83]), as their properties are related to those of the double-box integrals considered here and in refs. [79,80].
The generalized discontinuity operators whose contours are given by eqs. (3.16) and (3.22), along with similar results from ref. [80], can be applied directly to computations of two-loop amplitudes in both numerical and analytic forms. Their construction also hints at deeper connections to the algebraic geometry of the corresponding Feynman integrals. where