Uniqueness of two-loop master contours

Generalized-unitarity calculations of two-loop amplitudes are performed by expanding the amplitude in a basis of master integrals and then determining the coefficients by taking a number of generalized cuts. In this paper, we present a complete classification of the solutions to the maximal cut of integrals with the double-box topology. The ideas presented here are expected to be relevant for all two-loop topologies as well. We find that these maximal-cut solutions are naturally associated with Riemann surfaces whose topology is determined by the number of states at the vertices of the double-box graph. In the case of four massless external momenta we find that, once the geometry of these Riemann surfaces is properly understood, there are uniquely defined master contours producing the coefficients of the double-box integrals in the basis decomposition of the two-loop amplitude. This is in perfect analogy with the situation in one-loop generalized unitarity. In addition, we point out that the chiral integrals recently introduced by Arkani-Hamed et al. can be used as master integrals for the double-box contributions to the two-loop amplitudes in any gauge theory. The infrared finiteness of these integrals allow for their coefficients as well as their integrated expressions to be evaluated in strictly four dimensions, providing significant technical simplification. We evaluate these integrals at four points and obtain remarkably compact results.


Introduction
The study of scattering amplitude in gauge theories is a fascinating subject, partly because of the close link it provides between theory and experiment, and partly because it continues to unravel interesting structures in quantum field theory. Our understanding of computing amplitudes has undergone a revolution over the past decade and a half, owing in large part to the development of on-shell recursion relations for tree-level amplitudes [1,2] and a purely on-shell formalism for loop-level amplitudes, the modern unitarity method . These powerful modern methods have to a large extent replaced the more traditional Feynman diagrammatic approach for tree-level and one-loop amplitudes. Thus, the current frontier is the development of systematic approaches for computing two-loop amplitudes. Impressive calculations of two-loop amplitudes have been done by means of Feynman diagrams, including all parton-level amplitudes required for e + e − annihilation into three jets, to give just one example [36,37] (these computations were subsequently used to extract α s to 1% accuracy from the three-jet LEP data [38,39]). The focus of this paper is to explore a different approach, the unitarity method at two loops. This method has proven very successful at one loop where it has rendered a number of amplitude calculations possible, in particular of processes with many partons in the final state. In this formalism, the one-loop amplitude is written as a sum over a set of basis integrals, with coefficients that are rational in external spinors, Amplitude = j∈Basis coefficient j Integral j + Rational . (1.1) The process-dependence thus resides in the integral coefficients which are the object of calculation within the unitarity-based approach. The determination of the coefficients is done by applying to both sides of this basis decomposition a number of cuts, defined in the basic variant of unitarity as computing the branch cut discontinuities across the various kinematic channels. As a result, the left hand side of eq. (1.1) is, by the Cutkosky rules, turned into a product of tree-level amplitudes, enabling the computation of one-loop amplitudes from tree-level data.
One-loop unitarity also exists in a more recent version, called generalized unitarity, in which the operation of taking cuts does not have any known interpretation in terms of branch cut discontinuities. Rather, generalized cuts are defined as a change of the integration range away from the real slice R D (where D = 4 − 2 ) into a contour of real dimension 4, embedded in C 4 . The resulting contour integrations compute residues that are unique to each of the basis integrals in eq. (1.1), enabling a direct extraction of their coefficients, using only tree-level amplitudes as input data.
Unitarity has also been applied beyond one loop, taking as the starting point a decomposition, similar to eq. (1.1), of the desired amplitude into a (typically overcomplete) basis and requiring agreement between the two sides on all cuts to determine the integral coefficients on the right hand side. Several impressive calculations have been done in this way, primarily in N = 4 supersymmetric Yang-Mills theory and in N = 8 supergravity. However, calculations of this nature require crafty choices of bases, and it is fair to say that no systematized use of generalized unitarity exists beyond one loop.
In ref. [34], the first steps were taken in developing a fully systematic version of generalized unitarity at two loops. In the approach followed there, the two-loop amplitude is decomposed as a linear combination of basis integrals, in similarity with eq. (1.1). The integral coefficients are determined by applying to both sides of this two-loop basis decomposition so-called augmented heptacuts, defined as a change of the integration range away from the real slice R D × R D into contours of real dimension 8, embedded in C 8 . These contours are particular linear combinations of the tori encircling the leading singularities of the integrand, and satisfying the consistency condition that any function which integrates to zero on the real slice R D × R D (where D = 4 − 2 ) must also integrate to zero on the C 8 -embedded contour. This constraint on the contour ensures that two Feynman integrals which are equal, possibly through some non-trivial relations, will also have identical maximal cuts. As explained in ref. [34], contours satisfying this consistency condition, or master contours, are guaranteed to produce correct results for scattering amplitudes in any gauge theory. A closely related, but distinct, approach is that of refs. [40,41], in which the heptacut integrand is reconstructed by polynomial matching in similarity with the OPP approach [21].
A perplexing feature of the contours obtained in ref. [34] is that they were not found to be unique, in contradistinction with the situation found at one loop [9,23]. In this paper, we will extend the results of ref. [34] to an arbitrary number of external legs, and at the same time show that the contours are actually unique, once proper identifications of the leading singularity cycles are taken into account. In this way, the situation at two loops becomes entirely analogous to the situation at one loop. This paper is organized as follows. In section 2, we introduce notation and formulas used throughout the paper. In section 3, we give our classification of the kinematical solutions to the maximal cut constraints of the general double-box integral and discuss the singularities of the Jacobian arising from linearizing these constraints. In section 4, we explain how to rephrase the problem of solving on-shell constraints as a geometric problem in momentum twistor space. In section 5, we prove that the master contours extracting the double-box coefficients of two-loop amplitudes are uniquely defined, once the sharing of Jacobian poles between kinematical solutions is properly taken into account. We then show that the doublebox-topology basis elements can be chosen to have chiral numerator insertions. In section 6, we give a detailed derivation of the analytic expressions of these chiral double-box integrals. In section 7, we provide our conclusions and suggest directions for future investigation.

Maximal cut of the general double box
Unitarity-based computations of two-loop amplitudes take as their starting point the expansion of the amplitude into a basis of linearly independent two-loop integrals, (2.1) The form of the right hand side is obtained by applying integral reductions to the Feynmandiagrammatic expansion of the amplitude, and the integrals Int i are referred to as master integrals (the rational contributions will disregarded in this paper). The process dependence resides in the integral coefficients c i , and the goal of unitarity calculations is to determine these coefficients as functions of the external momenta.
In this paper, we will be concerned with the coefficients of the integrals in eq. (2.1) containing the maximal number of propagators. These integrals turn out to have the doublebox topology, 1 illustrated in figure 1. They are defined by where D = 4 − 2 . The function Φ( 1 , 2 ) multiplied into the integrand will be referred to as a numerator insertion. The special case of Φ = 1 produces a so-called scalar double-box integral. The right hand side of eq. (2.1) will, depending on the number of external momenta in the process in question, contain several integrals of the double-box topology with various numerator insertions (for four external momenta, there are two master integrals; the two used in ref. [34] involved the insertions Φ = 1 and Φ = 1 · k 4 ). In order to determine the integral coefficients in eq. (2.1), one applies to both sides a number of so-called generalized-unitarity cuts. These can roughly speaking be understood as a replacement of the seven propagators in eq. (2.2) by seven δ-functions whose arguments are given as the corresponding inverse propagators. These δ-functions thereby solve what are called the on-shell constraints, whose solutions 1 , 2 are generically complex. The effect of applying such generalized cuts to eq. (2.1) is to turn the loop integrals into contour integrals in the complex plane. Then (roughly speaking) by choosing the integration contours to encircle poles unique to each master integral in this basis decomposition, one may extract their coefficients, thereby determining the amplitude. We call such contours master contours and will discuss them in much greater detail in section 5.
To solve the on-shell constraints (2.3)-(2.9) for an arbitary number of external momenta, it 1 When the number of external states exceeds four, the leading topology is that of a pentagon-box or a double-pentagon. However, we expect the coefficients of such integrals to be simpler to extract due to the explicit octacuts they contain.
proves useful to introduce, following refs. [21,23], null vectors K µ 1 and K µ 2 that lie within the plane spanned by K µ 1 and K µ 2 . Note that the vectors K µ 1 and K µ 2 are not necessarily assumed to be null, but the vectors K µ 1 and K µ 2 are appropriate null linear combinations. Similarly, we introduce vectors K 4 and K 5 in the plane spanned by K 4 and K 5 . Using these vectors, a convenient parametrization of the loop momenta is given by Re-expressed in terms of the loop momentum parametrization (2.10)-(2.11), the on-shell constraints (2.3)-(2.8) (corresponding to cutting the six outer propagators in figure 1) take the form and γ 2 is defined in analogy with γ 1 . We refer to refs. [21,23] for more details. We see that in this parametrization, the variables α 1 , α 2 , β 1 , β 2 are directly fixed, while the remaining variables obey simple constraints of the form α 3 α 4 = constant. On the solution of the above equations, cutting the central propagator gives rise to the equation Cutting the seven propagators visible in the double-box graph in figure 1 will only fix seven out of the eight integration variables in the two-loop integration. Nevertheless, as was explained in ref. [42], due to the Jacobian factors arising from solving the δ-functions, the measure for the remaining variable can develop poles at specific locations. The last integration can then be performed on a small circle enclosing such poles, effectively pulling out an eighth propagator. Such poles are referred to as the leading singularities of the integrand.
For future reference, let us provide a few details of the computation of the maximal cut of the double-box integral in eq. (2.2). Replacing all propagators in eq. (2.2) by δ-functions and integrating out the six that correspond to the outer propagators in figure 1, the heptacut measure on a given kinematical solution is . (2.14) We note that the variables α 3 or β 3 are not always good integration variables: on certain solutions to eq. (2.12) they may happen to be constant. In such cases, they should be traded for α 4 or β 4 through dα 3 α 3 → − dα 4 α 4 , and/or dβ 3 β 3 → − dβ 4 β 4 . Notice the relative signs, which are essential to ensure the global consistency of the residues. These arise from the fact that, in solving for the δ-functions, one should use determinants, not absolute values of determinants (see, for instance, the discussion in ref. [43]).
Integrating out the remaining δ-function in eq. (2.14), for each kinematical solution, produces the corresponding maximal cut where z ≡ α 3 and (2.17) Here α 1 , α 2 , α 4 (z), β 1 , β 2 are given by eq. (2.12), whereas z is unconstrained, reflecting the degree of freedom left over after imposing the seven cut constraints. Similar formulas arise when solving instead for z = α 4 , β 3 or β 4 . Despite appearances, we will see that in all cases with less than 10 massless particles, the argument of the square root in eq. (2.15) is in fact a perfect square.
Our goal in the next section will be to determine when and where does the integrand of eq. (2.15), referred to throughout this paper as the (heptacut) Jacobian, give rise to poles. As we will find, these poles are naturally associated with three-point vertices in the double-box graph illustrated in figure 1, and their locations can be understood in a simple way.

Kinematical solutions and Jacobian poles
In this section we consider the classes of solutions to the joint heptacut constraints (2.3)-(2.9). As the loop momenta have a total of eight degrees of freedom (α 1 , . . . , α 4 , β 1 , . . . , β 4 ), the result of imposing seven on-shell constraints will be to fix all but one of these parameters. The various choices of freezing the loop parameters to particular values that solve these constraints span a number of distinct kinematical solutions, whose unconstrained variable z ∈ C parametrizes a Riemann surface (for example, a Riemann sphere). As we shall see, the Riemann surfaces associated with the kinematical solutions are not disjoint, but rather they have pointwise intersections located at the poles of the Jacobian discussed in the previous section.
The number of kinematical solutions to the heptacut constraints is determined by the distribution of external momenta at the vertices of the double-box graph, and an important role in the classification is played by the vertices that join three massless lines. In order Figure 1. The general double-box integral. The · · · dots at each vertex represent the presence of an arbitrary number of massless legs. Each of the vertices, shown as gray blobs, is given a label i = 1, . . . , 6 which equals the index of the associated external momentum K i .
to state the classification, we introduce some notation which will be used throughout the paper, for j = 1 for j = 2 for j = 3 . (3.1) The variable µ j keeps track of whether each of the three vertical lines in the double-box graph in figure 1 is part of some three-point vertex or not, and respectively equals zero or one. For mnemonic convenience, we will denote the values of µ j by letters as follows Finally, the notation I (N 1 ,N 2 ,N 3 ,N 4 ,N 5 ,N 6 ) will be used to refer to a double-box integral with N i external massless legs attached to vertex i.
To give an example of how three-point vertices play a role in determining the number of kinematical solutions, let us consider the third equation in eq. (2.12), We observe that if the right hand side is nonzero, one gets an invertible relation between α 3 and α 4 , leaving either of them as an equivalent free parameter. If, on the other hand, the right hand side is zero, this equation has two distinct solutions, α 3 = 0 or α 4 = 0. The latter situation occurs whenever the leftmost vertical line of the double-box graph is part of some three-point vertex (in the above notation denoted by m), and the splitting of one into two solutions of eq. (2.12) is a reflection of the existence of two types of massless on-shell three-point vertices in 3 + 1 dimensions [44]. Indeed, assuming for the moment that S 1 = 0, it follows from eq. (2.10) that α 3 = 0 implies Thus we see that the three-point vertex is special in that it connects momenta with the property that either, as in eq. (3.4), their square-bracket spinors are aligned, or, as in eq. (3.5), their angle-bracket spinors are aligned. In the following, we shall respectively denote the two cases in eqs. (3.4) and (3.5) with a and ⊕ label, and refer to a label as the chirality of the vertex. Unlike a three-point vertex, a vertex with four or more legs does not have a well-defined chirality.
Below we discuss the number of kinematical solutions to the heptacut constraints (2.3)-(2.9) and the intersections of their associated Riemann surfaces for each of a total of four cases. These four cases are defined by having all, exactly two, exactly one and none of the vertical lines in the double-box graph be part of some three-point vertex. Let us first consider the integral topologies where each vertical line of the double-box graph is part of at least one three-point vertex. This category includes, for instance, all topologies with four or five massless external states, but also an infinite sequence of topologies at higher points.
The case of four external massless states was studied in detail in ref. [34] where it was shown that the number of kinematical solutions to the heptacut constraints (2.3)-(2.9) is six -as explained there, the solutions are uniquely characterized by the distribution of chiralities (⊕ or ) at the three-point vertices. Disallowed distributions are those with uninterrupted chains of same-chirality vertices along the vertical or horizontal lines; all other configurations are allowed for generic external kinematics.

Classification of kinematical solutions
As it turns out, the classification of kinematical solutions in the case of four massless external states extends uniformly to cover all case 1 topologies. We can establish this in two steps.
Let us start by noting that for an allowed solution, three-point vertices lying on a common vertical line must have opposite chirality. For the leftmost and rightmost lines, this is already visible from the discussion around eqs. (3.4) and (3.5). For example, having two ⊕ vertices on the leftmost vertical line cannot be achieved for generic external momenta because it would require and thus K 1 ·K 2 = 0. Similarly, for generic external momenta, forcing two ⊕ vertices to appear in the central vertical line can be shown to require the four-momentum in the central propagator to vanish, leaving no free parameter.
This already places an upper bound of 2 3 = 8 on the number of kinematical solutions. Not all of these configurations are allowed, however. For instance, three ⊕ vertices lying on a horizontal line would force all angle bracket spinors on that line to be proportional to each other, in analogy with eq. (3.6).
More generally, consider the distribution of chiralities at the vertices of the 7-particle topology shown in figure 2. This distribution of chiralities does not obviously impose any constraints on the external kinematics. Closer inspection, however, reveals that this assignment contains a one-loop sub-box on the right loop (of the two-mass-easy type), with opposite chiralities at its two massless corners. But as is known from studies of one-loop boxes, for generic external momenta, opposite corners of a two-mass-easy box must have identical chiralities [9]. Thus, the configuration in figure 2 in fact does impose constraints on the external momenta and is therefore disallowed. We have verified, by exhaustion, that all chirality assignments not forbidden in such ways lead to healthy kinematical solutions.
To summarize this discussion, we have derived three simple rules which establish that there are exactly six kinematical solutions for all topologies within case 1: • Rule 1. Two same-chirality vertices cannot appear in a vertical line.
• Rule 2. Three same-chirality vertices cannot appear in an horizontal line.
• Rule 3. In rule 2, opposite-chirality vertices at opposite corners of a one-loop sub-box should be counted like same-chirality adjacent vertices (cf. the right loop in figure 2).

Interpretation of Jacobian poles
As it turns out, the six kinematical solutions to the heptacut constraints (2.3)-(2.9) are not completely disjoint: as illustrated in figure 3, for any given kinematical solution S i , depicted there as a Riemann sphere, there are two special points where it coincides with a different kinematical solution. We now proceed to locate these special points and show that the six Riemann spheres link into a chain. To clarify the exposition, we will consider a particular representative of case 1 and write out explicitly the kinematical solutions and the Jacobian determinants. For example, let us choose the integral topology whose vertex momenta (as defined in figure 1) are with the k i being lightlike vectors. In terms of the spinor ratios the six kinematical solutions S 1 , . . . , S 6 are obtained by fixing the parameters of the loop momenta to the values S 1 , . . . , S 6 : α 1 = 1 , β 1 = 0 α 2 = 0 , β 2 = 1 (3.10) S 1 : with z ∈ C a free parameter. The associated heptacut Jacobians are (3.12) In the first two lines of eq. (3.12), the plus or minus signs refer respectively to the first or second indicated kinematical solution. Now, consider first the intersection between S 4 and S 6 . In S 4 we have α 3 = 0 but α 4 free, while in S 6 we have α 4 = 0 and α 3 free; the intersection is simply a point, located in S 4 at α 4 = 0 and in S 6 at α 3 = 0. At this point, β 3 equals zero while β 4 takes on a finite value explicitly given in eqs. (3.9) and (3.11).
To understand better why S 4 and S 6 coincide at a point, let us examine what is happening to the loop momentum 1 at z = 0 in S 4 . By assumption, either vertex 1 or 2 in figure 1 is a three-point vertex; let us consider here the former case. It is straightforward to see that in the parametrization (2.10) the on-shell constraints (2.3)-(2.5) and (2.9) are solved within S 6 by setting (α 1 , α 2 , α 3 , α 4 ) = ( S 2 +γ 1 γ 1 , 0, z, 0). At z = 0 we then observe that i.e., the loop momentum is collinear with that of a massless external particle. This collinearity can be immediately understood from figure 3: at the intersection of S 4 and S 6 , the lower-left vertex must simultaneously be of the and ⊕ type, and therefore the momenta connected by this vertex are mutually collinear. Moreover, when both of the momenta K 1 and K 2 are massless, the simultaneous collinearity conditions at the two leftmost vertices imply that the momentum of the particle exchanged between these vertices must vanish. Indeed, we see that in this case, eq. (3.13) implies that 1 −K 1 = 0. Physically, this corresponds to a soft divergence region, giving rise to an infrared divergence in the original two-loop integral. In a gauge theory, the exchanged soft particle producing such a singularity will necessarily be a soft gluon, as can be argued from the behavior of the three-point vertices.
Similarly, the intersection between S 2 and S 5 occurs at a point where α 3 = α 4 = 0. By symmetry, there are similar intersections at points where β 3 = β 4 = 0, merging S 1 with S 6 , and S 3 with S 5 . Finally, the intersections between S 1 and S 2 , and between S 3 and S 4 , occur at points where the momenta in the central three-point vertices become collinear.
Let us conclude this discussion by the observation that the poles of the Jacobian determinants in eq. (3.12) coincide with the intersection points of the six Riemann spheres shown in figure 3. We have checked that this phenomenon extends to all integral topologies in case 1.

Residue relations across solutions
At the location of any Jacobian pole, each loop momentum i as evaluated from either of two intersecting spheres assumes identical values; for example, As a result of this, given an arbitrary function f ( 1 (z), 2 (z)) that does not share these poles, one has the identities Res z=Q • singularities in refs. [42,45,46]. Pleasingly, we find that they are directly related to the physical collinear and infrared singularities of the theory. This is the case where exactly two vertical lines of the double-box graph are part of some three-point vertex. The analysis leading to rules 1-3 in section 3.1.1 remains valid and shows that there are exactly four kinematical solutions for any topology belonging to this case. These solutions are uniquely characterized by those assignments of vertex chiralities where no two ⊕ or occur in the same vertical line. Based on our insights from the previous subsection, it is natural to expect that the four Riemann spheres associated with the kinematical solutions again have pointwise intersections and are linked into a chain. This picture turns out to be correct, and we will now elaborate on some of its details.
For the double-box topologies of the type (m, M, m), one finds that the intersection points coincide with the poles of the Jacobian in eq. (2.15), in complete analogy with case 1. Moreover, each of these intersection points is associated with the simultaneous collinearity of the momenta in some three-point vertex of the double-box graph, again in exact analogy with case 1. The kinematical solutions of subcase (m, M, m) are illustrated in figure 4(a).
Seemingly, a new technical issue arises for the topologies of type (M, m, m): 3 the square root in eq. (2.15) suggests that the Jacobian contains branch cuts. Despite appearances, the radicand is in fact a perfect square, as can easily be seen by writing the central-propagator Plugging eq. (3.21) into eq. (2.14) and performing the last integration yields an explicitly rational formula for the Jacobian, containing poles rather than branch cuts. The (M, m, m) subcase thus presents no new features compared to the (m, M, m) subcase, except in one regard, illustrated in figure 4(b): the number of points in each Riemann sphere at which one of the loop momenta becomes infinite. This can be understood as follows. For solutions S 2 and S 4 , the fact that S 1 S 2 = 0 implies that α 3 ∼ z and α 4 ∼ 1 z , allowing for two distinct points on each of these spheres at which the left loop momentum 1 becomes infinite, denoted respectively as ∞ L,1 and ∞ L,2 in figure 4(b). On the other hand, in solutions S 1 and S 3 the loop momentum 1 assumes a constant value independent of z. 4 In particular, neither of these spheres contain any point at which 1 becomes infinite.
The above reasoning readily extends to case 1 and explains the positioning of the infinity poles in figure 3.
From figures 3 and 4 we observe that in both case 1 and 2, the number of independent residues one can take is 8. Here, the qualifier "independent" refers to the fact that on any given Riemann sphere, the residues necessarily add to zero, allowing any one residue to be expressed in terms of the remaining ones on the sphere. Thus, counting the number of poles shown in figures 3 and 4, and subtracting the number of Riemann spheres to compensate for the redundancy, we find 8 independent residues in all cases considered so far. In double-box topologies where exactly one vertical line of the graph is part of some threepoint vertex, the rules of section 3.1.1 imply that there are two kinematical solutions.
In the (M, M, m) case, one of the two Riemann spheres can be parametrized by z ≡ α 3 (and its parity conjugate by z ≡ α 4 ), and so the equations (2.15)-(2.18) readily apply. Because B −1 = 0, the Jacobian (2.15) is manifestly a rational function and has only poles 3 By assumption, at least one of the vertex momenta K3 or K6 (defined in figure 1) vanishes. Without loss of generality, we take K6 to vanish here. 4 To be somewhat more detailed, the chirality distributions in solutions S1 and S3 allow us to construct an "effective" one-loop box, obtained by collapsing the horizontal propagators between same-chirality labels on the right loop. As it turns out, the solution to the quadruple-cut constraints of this one-loop box is exactly equal to the left loop momentum 1 in the original heptacut double box. But as the quadruple cut freezes all components of the one-loop box momentum, this implies that the 1 obtained from S1 and S3 cannot have any dependence on z. in z. Exactly as in the previous cases, these poles are located at the intersections of the Riemann spheres; in particular, the Jacobian has exactly two poles on each sphere. In the (M, m, M) case, assuming again (without loss of generality) that K 6 = 0, we can proceed as in the (M, m, m) case above, following eq. (3.21). The same expression remains valid here and makes manifest the fact that the Jacobians are rational functions (of the variables α 3 , α 4 , β 3 or β 4 ). Again, the poles of the Jacobians coincide with the intersections of the Riemann spheres corresponding to the kinematical solutions.
The This is the case in which the double-box graph contains no three-point vertices. For the scattering of massless particles, the first time this occurs is for 10 particles, as depicted in figure 6.
To analyze this case, we return to the Jacobian determinant in eq. (2.15) which takes the form is a quartic polynomial. Numerically, we find that for generic 10-particle kinematics, the four roots r i of this polynomial are distinct. This means that, contrary to the previous cases, the Jacobian contains genuine branch cuts (meaning that they cannot be removed by any redefinition of z). The integration variable z in eq. (3.22) therefore parametrizes a two-sheeted cover of the Riemann sphere. This is topologically equivalent to an elliptic curve (i.e., a genus one Riemann surface), as Figure 6. The integral I (2,2,1,2,2,1) , the simplest example of an integral with more than three particles at all vertices, and whose heptacut Jacobian J(z) accordingly has branch cuts that cannot be removed by any reparametrization z → ϕ(z). As argued in the main text, this is presumably related to the appearance of functions in the analytic expression for I (2,2,1,2,2,1) which cannot be expressed in terms of generalized polylogarithms.
illustrated in figure 7. In particular, there is a single class of kinematical solutions to the heptacut constraints (2.3)-(2.9) in this case. For an elliptic curve there are two natural cycles over which the z integration in eq. (3.22) can be performed, generalizing the notion of a residue -namely, its topological cycles Γ 1 and Γ 2 , respectively shown in red and blue in figure 7. In terms of the z variable, these are cycles which enclose a pair of branch points. Integrations over such cycles produce so-called complete elliptic integrals of the first kind K(t) where the argument t is some cross-ratio of the four roots of the radicand Q(z). As they arise when performing the loop integration on a compact T 8 contour, the integration cycles Γ 1 and Γ 2 define leading singularity cycles of the double-box integral.
As illustrated in figure 7, the number of poles at which one of the loop momenta becomes infinite is 8. This can easily be explained as follows. The fact that S 1 S 2 = 0 implies that α 3 ∼ z and α 4 ∼ 1 z , allowing for two distinct points on each sheet of the elliptic curve at which the left loop momentum 1 becomes infinite; these points are denoted as ∞ L,i in the figure. Moreover, since each of the sheets can equivalently be parametrized in terms of β 3 or β 4 , the fact that S 4 S 5 = 0 implies that β 3 ∼ z and β 4 ∼ 1 z , allowing for two distinct points on each sheet at which the right loop momentum 2 becomes infinite; these points are denoted as ∞ R,i . Thus, there are in total 8 poles on the elliptic curve.
The residues of these 8 poles are not all independent, however. For instance, their sum is zero since it corresponds to a contractible cycle, as can be seen from figure 7. If all the infinity poles were only simple poles, which would be the case if all numerator insertions in the double-box integral were at most linear in each of the loop momenta 1 , 2 , there would exist a second, less obvious relation 5 . However, the large-momentum behavior of theories 5 This relation is easier to describe when the elliptic curve is viewed as the complex plane modulo the doubly-periodic identification z z + 1, z z + τ . While the first relation mentioned in the main text arises from integrating a form ω(z) along the boundary of a fundamental domain, the second relation arises Here, the Jacobian develops branch cuts, and the heptacut loop-momentum parameter z thus parametrizes a two-sheeted cover of the Riemann sphere. The two sheets are shown glued together along their branch cuts, illustrating how this Riemann surface is topologically equivalent to an elliptic curve. The topological cycles Γ 1 and Γ 2 , respectively shown in red and blue, enclose two distinct pairs of branch points and provide natural contours of integration for the heptacut double box in eq. (3.22). We observe that there are eight poles at which one of the loop momenta becomes infinite. In this case one finds a total of nine independent leading-singularity cycles, as explained in the main text.
such as pure Yang-Mills or QCD are not such as to produce only simple poles, and so this relation does not apply. Including the two topological cycles Γ i in our counting, we thus find a total of 9 independent leading-singularity cycles in case 4, in contradistinction with the previous cases 1-3.
Finally, let us summarize our discussion of the number of solutions to the maximal cut of the double-box graph in figure 1. We have observed that as the number of three-point vertices in the double-box graph grows, the number of associated Riemann surfaces increases from one surface (of genus one) to two spheres linked by pointwise intersections, and further on to four, and finally six spheres thus linked. The branching of kinematical solutions can be understood intuitively from the observation that the equation x 1 x 2 = m with m = 0 has a single connected component as its solution when x 1 and x 2 are allowed to be complex; but two essentially-disconnected components when m = 0. Applied to, e.g., the equation α 3 α 4 ∝ S 1 S 2 in eq. (2.12), this insight leads us to expect a splitting of one Riemann sphere into two as S 1 S 2 → 0. 6 This is indeed what happens, as exemplified by the splitting of the sphere S 1 in figure 5 into the spheres S 3 and S 4 in figure 4(a) whose left-most pair of double-box vertices then acquires chiralities. Taking the limit µ 2 → m of figure 4(a), the middle pair of vertices will acquire chiralities, splitting the spheres S 2 and S 4 into the respective pairs (S 1 , S 2 ) and (S 3 , S 4 ) in figure 3. In contrast, the solutions (S 1 , S 3 ) from integrating ω(z)z. In the absence of double poles, the latter relation relates the sum of the residues weighted by z to integrals over Γ1 and Γ2. We refer the reader to ref. [47] for more details; in particular, to Chapter III, Proposition 2.1. 6 We thank D. Kosower for this observation.
in figure 4(a) admit only one chirality assignment to the middle vertices, as dictated by the rules in section 3.1.1, and are transformed into the solutions (S 5 , S 6 ) in figure 3. As the resulting six solutions shown in figure 3 have chiralities assigned to all pairs of vertices, the splitting of Riemann spheres terminates at this stage. Let us finally remark that by giving generic small masses to the internal lines of the double-box integral with four lightlike external momenta, we expect that the six spheres that arose in the massless case are turned into a smooth elliptic curve. It would be interesting to study the situation for general patterns of internal masses.

Maximal cuts versus integrated expressions
This subsection lies outside the main scope of this paper and could be omitted on a first reading.
Widely propagated folklore holds that there should be a close connection between the maximal cut of a given integral and the analytic form of its integrated expression. The appearance of an elliptic curve in case 4 provides closer evidence of such a connection, as we will now argue. As shown in ref. [48], eq. (8.1), the integral in figure 6 can be represented as a one-scale integral as follows whereQ(u ) defines the same elliptic curve as Q(z) in eq. (3.22). It may be shown that the sunrise integral with massive propagators admits a very similar integral representation, with the integrand of eq. (3.23) containing log(· · · ) instead of Li 3 (· · · ) [49]. This integral was studied analytically in great detail in refs. [50,51] and references therein, and was found not to be expressible in terms of polylogarithms. Given the similarity with eq. (3.23), we thus find it extremely unlikely that I (2,2,1,2,2,1) is expressible in terms of (multiple) polylogarithms. Thus the topology of the Riemann surface parametrizing the heptacut solutions appears to be reflected in the integrated expression. Moreover, as we argue in appendix B, such more general functions must necessarily be present in the scattering amplitudes of N = 4 super Yang-Mills theory. In this appendix, we point out that a particular two-loop N 3 MHV amplitude for the scattering of 10 massless scalars is given precisely by the integral I (2,2,1,2,2,1) alone. This suggests that the realm of N = 4 SYM extends beyond that of polylogarithms.

Twistor geometry of two-loop maximal cuts
Momentum twistor space provides an appealing geometric rephrasing of the problem of setting massless propagators on-shell. Accordingly, we devote this section to explain the momentum twistor geometry of two-loop maximal cuts, the one-loop case having already been presented in ref. [52]. Besides geometric elegance, this formulation has the conceptual advantage of being inherently coordinate free. However, as the results of this section are not used elsewhere in this paper, this section can be skipped in a first reading.
The main notions are briefly reviewed in section 4.1, after which we describe the twistor space geometry of the heptacut in section 4.2.

Generalities
Given a sequence of null momenta k 1 , . . . , k n , where k µ i = 1 2 k i |γ µ |k i ], momentum twistors are four-component objects defined as [53] The definition is invertible: given a configuration of momentum twistors, a simple and explicitly known formula recovers the four-momenta k i [53]. An often useful fact is that any sequence of momentum twistors gives on-shell momenta which respect momentum conservation. In other words, the momentum twistors are free variables which solve all phase-space constraints. However, in the following we will not need this fact, as we will take the momenta k i to be given and construct momentum twistors out of them.
In analogy with the spinors k i |, the momentum twistors Z a i are defined only up to an overall rescaling; that is, they are points in three-dimensional complex projective space CP 3 . Due to this projective invariance, the subspace spanned by two momentum twistors Z i−1 and Z i defines a line in momentum twistor space which we will denote (i−1 i). This two-dimensional span is naturally recorded by the 4 × 2 matrix (Z i−1 Z i ). But since only the two-dimensional subspace itself matters, rather than any basis chosen for it, we can multiply this matrix from the right by an appropriate GL(2) matrix to obtain (assuming where σ αβ µ are the 2×2 Pauli matrices. Observe that the line in eq. (4.2) encodes the value of the associated region momentum x i , thereby identifying lines in momentum twistor space with points in region momentum space. Conversely, given a region momentum, one can construct the matrix in eq. (4.2) and interpret it as a line in momentum twistor space.
The kinematic invariants that can arise when considering planar graphs take the form [53], whose verification we omit here, is that and abcd is the antisymmetric Levi-Civita symbol. What will be important in the following is the geometrical interpretation of eq. (4.3): geometrically, what eq. (4.3) implies is that (k i +k i+1 +· · ·+k j ) 2 equals zero if and only if the lines (i−1 i) and (j j+1) intersect in momentum twistor space CP 3 . Indeed, if i−1 i j j+1 equals zero, one of the four twistors Z i−1 , Z i , Z j , Z j+1 can be written as a linear combination of the others. This means that there is some point in CP 3 which is simultaneously a linear combination of Z i−1 and Z i on one hand, and of Z j and Z j+1 on the other. That is, the lines (i − 1 i) and (j j + 1) intersect. As noted above, each loop momentum gives rise to a line in CP 3 through its associated region momenta. Thus, when such lines intersect, propagators become on-shell.

Heptacut example
Let us see how this works for the 7-point topology shown in figure 2. For conciseness, we will only describe the kinematical solution in which the two three-point vertices in the lower row of the double-box graph have chirality and ⊕, respectively, and the chirality in the upper-right corner is ⊕ (corresponding to the kinematical solution S 5 , see figure 3). The other distributions of chiralities are treated in an entirely analogous way.
The first step is to draw the configuration of momentum twistors associated with the external data, as shown in figure 8(a). The rule is the following: each region in the exterior of the double-box graph gives rise to a line. Whenever two regions are separated by a massless external leg of the graph, the corresponding lines intersect.
The second step is to construct two lines (AB) and (CD), respectively corresponding to the left and right loop momentum, which have appropriate intersections with the lines encoding the external momenta. More precisely, imposing the heptacut constraints, the line (AB) must intersect with the three lines (71), (12) and (34); and the line (CD) must intersect with (45), (56) and (71). Furthermore, the lines (AB) and (CD) must intersect each other, corresponding to having the central propagator cut. Thus, the problem of setting propagators on-shell translates into a geometry problem involving straight lines in a threedimensional space.
Let us solve it in steps. First, (CD) has to intersect (45) and (56). There is one obvious solution: (CD) can pass through the point (5). Actually, there is another solution: (CD) could be any line inside the plane (456). These two discrete solutions correspond to the ⊕ and vertices, respectively. For example, the ⊕ vertex joins momenta whose holomorphic spinors k i | are proportional, implying equality of their corresponding momentum twistors (4.1). Thus, as the upper-right vertex of the double box is assumed to be ⊕, (CD) is a line passing through the point (5). Furthermore, (CD) must intersect (71). Since we have just described two points on the line (CD), and since only the line they span is meaningful, it is straightforward to write down a parametrization for (CD), where z is a free variable. Next, we consider (AB). As the lower-left vertex of the double box is assumed to be , (AB) must be inside the plane (712). This characterization of vertices follows from the above characterization of ⊕ vertices, given that parity interchanges points and planes in momentum twistor space [53]. Furthermore, (AB) must intersect with the line (34). This implies that (AB) must pass through the point where the line (34) intersects the plane (712). This information allows us to write down the following parametrization of (AB), (4.6) Indeed, the coefficients in eq. (4.6) are such that 712A = 712B = 0; that is, the line (AB) is contained in the plane (712). Finally, what remains to be achieved is to have the lines (AB) and (CD) intersect each other; this is realized by setting w = z. We have thus fully characterized the solution S 5 as a function of one free parameter z. Through eq. (4.2), one can proceed to extract the loop momenta 1 and 2 as functions of this variable.
Let us now consider the behavior of the solution as a function of z (see figure 8(b)). For instance, as z → 0, (B) approaches the point (1) Using only such geometrical reasoning, one may construct all the kinematical solutions to the heptacut constraints of a given double-box integral and proceed to prove that they link into a chain, as already discussed and proven in sections 3.1-3.4.

Uniqueness of two-loop master contours
In reference [34], the first steps were taken in extending maximal unitarity beyond one loop. It was argued there that the contours of integration associated with two-loop maximal cuts cannot be chosen arbitrarily; rather, they are subject to the consistency condition that any function which integrates to zero on the real slice R D × R D (where D = 4 − 2 ) must also integrate to zero on the maximal-cut contour. This constraint on the contour ensures that two Feynman integrals which are equal, possibly through some non-trivial relations, will also have identical maximal cuts. As argued in ref. [34], contours satisfying this consistency condition are guaranteed to produce correct results for scattering amplitudes in any quantum field theory. To be more accurate, the validity of the approach extends to all quantum field theories in which all states are massless.
In particular, no assumptions need to be made regarding the powers of loop momentum present in numerators. This is due to the nontrivial, but empirically true fact that the number of master integrals in a given topology is independent of the loop-momentum power counting beyond a certain threshold (see, for instance, ref. [54]). Therefore, this method can be applied indiscriminately to N = 4, 2, 1, 0 Yang-Mills theory, massless QCD, gravity or scalar theories (in the latter two cases, as soon as the analysis has been extended beyond the planar sector).
The maximal cut of the double-box integral is in general a contour integral in the complex plane. For the case of four massless external momenta k 1 , . . . , k 4 , it turns out that the contour may encircle 14 different leading singularities with some a priori undetermined winding numbers, in the notation of ref. where i = 1, . . . , 6, j = 5, 6 and χ ≡ s 14 s 12 . Here, the 12 winding numbers a 1,i and a 2,i are associated with Jacobian poles whereas the 2 winding numbers a 3,5 and a 3,6 are associated with the poles ∞ R in S 5 and S 6 in figure 3. The consistency conditions on the maximalcut contours were shown to translate into the following linear constraints on the winding numbers, a 1,2 + a 1,5 − a 1,4 − a 1,6 = 0 a 2,1 + a 2,2 − a 2,3 − a 2,4 = 0 a 2,6 − a 1,1 − a 2,5 + a 1,3 = 0 a 3,5 − a 3,6 = 0 a 1,2 + a 1,5 + a 1,4 + a 1,6 = −a 2,6 + a 1,1 − a 2,5 + a 1,3 + a 3,5 + a 3,6 a 3,5 + a 3,6 = − between tensor double-box integrals. Imposing these constraints on the winding numbers leaves 14 − 6 = 8 free parameters in the contours. On the other hand, in the case of four massless external momenta, there turn out to be exactly two linearly independent master integrals of the double-box topology. The heptacuts of the particular masters used in ref. [34] were found there to evaluate to where we remind the reader of the numerator insertion notation explained in eq. (2.2). As shown in detail in this reference, it is possible to find contours (a 1,i , a 2,i , a 3,j ) satisfying the constraint equations (5.2) with the additional property of setting the right hand side of eq. (5.4) to zero while setting the right hand side of eq. (5.3) to one (or vice versa). Such contours thus isolate the contribution of a single master integral in the basis decomposition (2.1) and are therefore referred to as master contours. With the 8 free parameters in the contours that remained in the previous paragraph, we thus find a total of 8 − 2 = 6 free parameters in the two-loop master contours. Extrapolating from the situation at one loop in which there is a unique contour associated with each of the basis integrals [9,23], we would not expect any free parameters in the master contours, and the appearance of 6 unconstrained variables comes as a surprise. By considering figure 3, this phenomenon can now easily be explained: all of the Jacobian poles belong to two Riemann spheres and so are counted twice in the above counting (5.1). Indeed, for example, the winding numbers a 1,2 and a 1,5 only arise in the combination (a 1,2 + a 1,5 ), as is visible in eq. (5.2). Due to this, a contour which encircles the z = 0 pole in S 2 and in S 5 with identical winding numbers in the two spheres, and no other poles, is equivalent to a zero-cycle. The addition of such zero-cycles defines an equivalence relation on the vector space spanned by the leading-singularity contours: any two contours related by the addition of such zero-cycles are equivalent. (a 2,6 , a 1,1 ) −→ (a 2,6 , a 1,1 ) + (ξ 6 , ξ 6 ) , (5.10) corresponding to the addition of a zero-cycle encircling each Jacobian pole with the winding numbers (ξ i , ±ξ i ) ∈ Z × Z on the two Riemann spheres containing the pole. These equivalence relations allow us to add to an arbitrary contour, characterized by the winding numbers (a 1,i , a 2,i , a 3,j ), the zero-cycle with (ξ i ) = (−a 2,1 , −a 1,2 , −a 2,5 , −a 2,3 , −a 1,4 , −a 2,6 ) to obtain an equivalent contour characterized by 8 independent parameters. This shows that the leading singularity contours are characterized by 8 rather than 14 winding numbers.

Invariant labeling of contours
To get around the redundancy built into the notation (5.1), we will from now on adopt the following notation for the independent winding numbers Here ω i∩j denotes the winding around the intersection point of S i and S j of a small circle supported in either of these spheres, with positive orientation in S i or with negative orientation in S j . Analogously, ω j,∞ R denotes the winding around the point ∞ R in S j where the right loop momentum 2 becomes infinite (see figure 3). Of course, we could trade some of the variables in eq. (5.11) for winding numbers around the other infinity poles in figure 3, but we find the above choice to be the most convenient. Also note that (5.11) has the added advantage over the notation (5.1) of not making reference to a particular parametrization of the Riemann spheres S i , hence the title of this subsection. The 8 winding numbers ω i in eq. (5.11) are equal to the following linear combinations of the a i,j Ω = (a 2,1 +a 2,2 , −a 1,2 −a 1,5 , a 2,5 −a 1,3 , a 2,3 +a 2,4 , −a 1,4 −a 1,6 , a 2,6 −a 1,1 , a 3,5 , a 3,6 ). (5.12) In analogy with the notation (5.11) for the winding numbers around the leading singularities, let us introduce the following notation for the residues. The residue at the intersection 8 Please note that in this section, in order to make the connection with ref. [34] as clear as possible, we adopt the conventions of this reference on the orientations of contours encircling poles in the various Riemann spheres and on the signs of the Jacobians. These conventions differ from those used in this paper, in particular, with respect to the minus signs described below eq. (2.14). The pattern of relative signs between the zero-cycle winding numbers (ξi, ±ξi) in eqs. (5.5)-(5.10) owes to the omission of these minus signs in ref. [34]. point of the spheres S i and S j , computed from the viewpoint of sphere S i , will be labeled Res i∩j . Alternatively, we could consider the same residue, but computed from the viewpoint of S j . As noted around eqs. (3.15)- (3.20), the result would be equal and opposite. Thus, we can re-express these identities in a very compact form by declaring Res i∩j to be antisymmetric, Res i∩j Φ(z) = −Res j∩i Φ(z) . Re-expressing the contour constraint equations (5.2) in terms of the winding numbers (5.11)-(5.12), they are found to take the form while the heptacut master double boxes used in ref. [34]  We observe that upon imposing the 6 constraint equations (5.14) on the 8 winding numbers ω i∩j given in eq. (5.11), we are left with 2 unconstrained parameters. This number of free parameters exactly equals the number of master double-box integrals at four points. In other words, we observe that in terms of the ω-variables, there is a unique master contour associated with each of the master double-box integrals in eqs. (5.15)-(5.16). These contours are respectively characterized by the winding numbers 0, 1, 1, 0, 1, 1, 1) and The observed uniqueness of master contours at two loops is in perfect analogy with the situation in one-loop generalized unitarity [9,23] and constitutes the main result of this section.

A basis with infrared finite master integrals?
As indicated in the basis decomposition of two-loop amplitudes in eq. (2.1), the integral coefficients are functions of the dimensional regulator . But in contrast to the situation at one loop, the O( ) contributions to the coefficients cannot be re-expressed as rational contributions to the amplitude, and these corrections therefore form an inevitable part of the two-loop integral coefficients. The physical significance of these O( ) contributions lies in the fact that in the basis expansion of the two-loop amplitude, they will multiply 1 k singularities in the integrated expressions for the two-loop integrals, thus producing finite contributions to the amplitude. Their extraction therefore poses an important problem. Unfortunately, as explained in ref. [55], the O( ) parts of two-loop integral coefficients are not obtainable from four-dimensional cuts. Instead, they must be computed by evaluating cuts in D = 4 − 2 dimensions, something which is technically much more involved. Ideally, one would like to circumvent the need for taking cuts in D = 4 − 2 dimensions, or at least limit such computations as much as possible. One way to achieve this would be to expand the two-loop amplitude in a basis that contains as many infrared finite integrals as possible. Indeed, although the expansion coefficients may still depend on , the physically relevant part of the coefficients multiplying IR finite integrals is purely O( 0 ) and may thus be obtained from strictly four-dimensional cuts. Of course, as two-loop amplitudes do have IR divergences of their own (as they necessarily must, to cancel the IR divergences of tree and one-loop amplitudes in the cross-section [56,57]), any basis of integrals must contain IR divergent integrals 9 . Nonetheless, it is plausible that by using a basis with a minimal number of IR divergent integrals one can minimize the work of extracting the physically relevant part of the basis integral coefficients. Focusing our attention to the double-box contributions to two-loop amplitudes, we thus turn to the question: can one find two linearly independent infrared finite integrals with the double-box topology?
A class of integrals with the property of infrared finiteness was introduced in ref. [58] where they were used to express the integrand of N = 4 super Yang-Mills amplitudes in a strikingly simple form. Here, we wish to investigate whether these so-called chiral numerator integrals can be used as master integrals for two-loop amplitudes in any gauge theory. For four lightlike external momenta there are, up to parity conjugation, two distinct chiral numerator integrals, which we can define as 10 where we remind the reader of the numerator insertion notation explained in eq. (2.2). Similar definitions can be given for an arbitrary number of external legs, replacing the spinors [i| or i| by their flattened counterparts [K i | or K i |. 9 Unless we are prepared to accept integral coefficients with 1 k singularities -which we are not. 10 We should stress that the chiral numerator integrals considered here are not exactly chiral in the sense of ref. [58], in which the general notion of a chiral integral was formulated in terms of the analytic structure of its integrand. There, chiral integrals were defined as integrals which have at most simple-pole singularities on each leading-singularity contour, and whose residues are all equal to either zero, or plus or minus one. In general, such integrals may differ from the chiral numerator integrals I+± considered here by the addition of integrals with fewer internal lines, such as triangle-boxes or double-triangles. Since we are not concerned with the latter integrals in this paper, this distinction will not be of relevance here, however.
The maximal cuts of these integrals are found to take the form I ++ cut = ω 6,∞ R (5.20) receiving contributions from a remarkably small number of leading singularities. The vanishing of a large number of leading-singularity residues reflects the absence of infrared singularities in the uncut integrals, as observed in section 3.1. Here, we could of course also have chosen to consider the integrals I −− and I −+ whose numerator insertions are obtained by parity conjugation ·| · |·] ←→ [·| · |· of eqs. (5.18)- (5.19). However, on the solution of the first four constraint equations of (5.14) (which express parity invariance of the contours), one finds that the maximal cuts of parity-conjugate integrals are equal. The largest potentially linearly independent set of chiral double boxes therefore consists, at four points, of two integrals which we choose as those given in eqs. (5.18)- (5.19). Considering these two integrals, let us now ask: are they linearly independent as master integrals? Equivalently, can one find two distinct contours, satisfying the constraint equations (5.14), each with the property of yielding a nonvanishing maximal cut for precisely one of these integrals?
Thus the integrals I ++ and I +− are linearly independent and may be used as master integrals for the double-box contributions to two-loop amplitudes in any gauge theory. Incidentally, as a bonus, the contours in eq. (5.22) take a somewhat simpler form than those in eq. (5.17) whose winding numbers display no easily discernible pattern.
The choice of using the chiral numerator integrals as master integrals provides a substantial simplification over other choices (such as that of eqs. (5.15)-(5.16)), as their infrared finiteness allows their expansion coefficients in the two-loop amplitude to be obtained from strictly four-dimensional cuts. But one might worry that the technical difficulty is simply shifted elsewhere, in particular to the analytical evaluation of these integrals. To counter this concern, we present in the next section a detailed analytical evaluation of these integrals and moreover find the result to take a remarkably compact form.

Analytical evaluation of chiral double boxes
In the previous section we found that the chiral numerator integrals in eqs. (5.18)-(5.19) form a basis of the master integrals with the double-box topology. The finiteness of these integrals allows for an economic extraction of their coefficients, which can be done directly in four space-time dimensions. In this section we turn to the question of evaluating these integrals analytically. We will focus our attention on obtaining analytical expressions in the special case of four lightlike external momenta, but we are hopeful that many of the ideas presented here will prove applicable for higher numbers of external legs as well.
Our procedure goes through several steps. First, we apply to the Feynman parametrized expression of the integrals a sequence of non-obvious changes of variables, motivated by recent work [48] on Mellin space transforms, but which make no reference to Mellin space. This enables us to obtain the symbols [59,60] of each of the chiral double boxes. The finiteness of the integrals, ensuring that they can be evaluated directly in four dimensions, is very helpful in this respect. Finally, imposing the constraints of uniform transcendentality, analyticity and Regge limits then allow us to unambiguously integrate these symbols, leaving us with the final expressions in eqs. (6.23)-(6.24).
Our starting point will be the (standard) Feynman parametrization formula, which gives for a double box with arbitrary numbers of legs 11 and x ij ≡ k i + · · · + k j−1 . The present notation is consistent with the dual coordinates x i introduced in section 4, setting x ij ≡ x i − x j . In eq. (6.1), we have dropped a number of terms of the form v 1 ·k 1 or v 1 ·k 2 , which vanish when v 1 , v 2 are chosen to correspond to the chiral numerators defined below eqs. (5.18)- (5.19).
In the special case of four lightlike external momenta, the second term on the last line of eq. (6.1) takes the form Furthermore, it turns out that This can be seen from eq. (6.1), because x 34 = 0 in this case (that is, K 3 = 0) and v 1 ·x 63 = 0, so only the b 5 term in the v 1 · i b i x i3 factor contributes. Due to similar simplifications in the v 2 factor, the first term is proportional to a 2 b 5 which is what the derivative produces. 11 We use the normalization I[· · · ] ≡ Plugging in the explicit expressions for v 1 and v 2 and including the appropriate constant of proportionality in eq. (6.4), one finds I ++ = −χ 2 1 + (1 + χ) ∂ ∂χ I 1 (χ) and The upshot is that we only have to compute I 1 . This particular simplification is probably specific to four points, but we welcome it.

Evaluation of I 1
It turns out to be possible to evaluate I 1 analytically, as we will now describe.

Step 1: Projectivization
We would prefer a form for I 1 in which such high powers of the integration variables do not occur within the denominator. This can be obtained using the following trick, which is not limited to the present integral. This might be called a "projectivization" trick. Thus, our goal in this subsection will be to derive the following equivalent form of eq. (6.3): The "1/vol(GL(1))" notation is an instruction to set any one of the variables a 1 , . . . , b I equal to 1 and integrate over the seven remaining ones; this is also why we wrote "d 7 " instead of "d 8 ". As the notation suggests, the result does not depend on which variable is set to 1, due to the scaling ("GL(1)") symmetry of the integrand. The reader not interested in following the derivation of eq. (6.6) from eq. (6.3) can safely skip to the next subsection.
The "trick" here is just a sequence of elementary, if non-obvious, changes of variables. First we use the δ-function in eq. (6.3) to perform for instance the a 1 integration, and make the change of variables 12 : (6.7) This brings eq. (6.3), after a convenient relabeling As the notation suggests, the result does not depend on which variable was used to solve the δ-function constraint. Second, we change integration variable c → c and then remove the prime, a 1 a 3 c + a 1 b 3 + a 3 b 1 + a 2 (6.9) Finally, we combine the three denominators in eq. (6.9) into one by introducing the two auxiliary Feynman parameters a I and b I ; these are such that integrating out a I and b I from eq. (6.6) produces eq. (6.9). (In a momentum twistor space computation, a I and b I would have appeared automatically as Feynman parameters for the propagators ABI and CDI , connecting to the so-called infinity point I. This is the origin of our notation a I , b I .) Our desire to obtain the form in eq. (6.6) was sparked by a recent paper [48] where very similar formulas for two-loop integrals were obtained with the help of Mellin space techniques. These formulas generically, as mentioned previously in eq. (3.23), involve an integral over a projective space with a quadratic form in the denominator, then integrated over one real variable (see, in particular, eq. (8.1) of ref. [48]). As emphasized by the authors of ref. [48], a great deal is known about such projective integrals, which makes these forms particularly attractive. As we will see below, such forms are indeed very convenient for computations. Here, our derivation used only elementary manipulations in Feynman parameter space 13 .

Step 2: Obtaining the symbol
Having obtained eq. (6.6), we observe that several of the integrations are trivial; for instance, we can immediately do the b 2 , b 3 and a 3 integrations, leaving We can easily integrate out one more variable, although this step will unavoidably produce a logarithm. For instance, doing the a 2 integration in eq. (6.10) yields a 1 c)) . (6.11) 13 We have found that the simple changes of variables described above, when applied to the double-box integral considered in eq. (8.1) of ref. [48], reproduces that formula exactly, giving an elementary derivation of it.
A pleasing surprise is that one can keep going like that, doing one integration at the time, and not hit any serious roadblock, besides the increasing length of the expressions. For instance, a 1 enters in eq. (6.11) at most linearly in all denominator factors and arguments of the logarithms, and so the a 1 integral can be done explicitly, producing dilogarithms.
If one continued like this, in the right order (doing b 1 and then c next), the integrals produced would take the same form as eq. (6.11), but with integrands involving transcendental factors of increasing transcendentality degree, for instance Li 3 , followed by Li 4 . This is manifest from the fact that each of these integrals take the schematic form × Li n (ratio of linear factors in t) + · · · (6.12) where the degree of the polylogarithm in the integrand is n = 2 or 3, depending on the integral. Integrations of this type raise the degree of transcendentality of the integrand exactly by one. It is relatively easy to check that the rational denominators involve only linear factors at all stages. Namely, as a general feature of such integrals, the rational part of the measure at a given stage is obtained simply by taking residues of the rational part of the measure at the previous stage. Finally, the a I integral (setting b I = 1 in this final stage) would be of the form (6.13) whose measure contains a squared linear factor. Using integration by parts this can be replaced by a boundary term plus an integral of the form da I (Li 3 (· · · )/rational + · · · ), both of which would manifestly give rise to polylogarithms of degree 4. Although the procedure outlined in the previous paragraphs is clearly feasible, carrying it out in practice would be cumbersome. It is much more efficient to perform the above integrals, exactly as described, but only at the level of the symbol. This can be done purely algebraically and will produce the symbol of I 1 . For instance, employing the algorithm described in appendix B of ref. [61], each step involves only solving linear equations and can be easily automated. The intermediate expressions are somewhat lengthy and will not be reproduced here, but in the end we obtain the very simple symbol, Knowing only this symbol, together with three other pieces of information, it is possible to reconstruct the function I 1 uniquely.

Step 3: Integrating the symbol
As we will now show, the symbol (6.14) can be integrated unambiguously by imposing the following three constraints on the integrated expression: 1. The fact that I 1 (χ) has the homogeneous transcendentality degree 4, as explained above. 14 2. The physical requirement that, on all physical sheets (which can be either −1 < χ < 0, −∞ < χ < −1 or 0 < χ < ∞, depending on the channel under consideration), the integral is analytic around χ = −1. This is because our integral, being planar, has a vanishing unitarity cut in the u-channel and hence, by the Cutkosky rules, cannot have a discontinuity in the u-channel.
The first constraint implies that where f 1 and f 2 are functions of homogeneous degree 4 whose symbol is consistent with eq. (6.14) and therefore must take the form Here "(symbol-free terms)" represent simpler transcendental functions multiplied by constants, such as π 2 Li 2 (· · · ) or ζ(3) log χ. The H's are harmonic polylogarithms [62]; explicitly, with 14 A technical issue which arises here is that, because of the integration by parts step, lower-degree contamination is a priori possible. This could arise if the rational part of the measure after the integration by parts step still contains a squared denominator 1/(aI + 1) 2 ; this would lower the transcendentality. We have verified at the level of the symbol of the (6.13) integrand that the integration by parts does not produce any such squared denominator. This means that such a squared denominator could arise only from a beyond-the-symbol ambiguity in the (6.13) integrand, hence would have to be explicitly proportional to π 2 or ζ(3). However, our Regge asymptotics constraints turn out to be strong enough to rule out such terms.
We now impose constraint 2, regularity at χ = −1. This has to be imposed separately on f 1 and f 2 because these transcendental functions multiply different rational prefactors. Because a function is analytic at a point if and only if its derivative is, the simplest way to proceed is to take a derivative. For instance, near χ = −1, we have and similarly for f 2 . This decomposition is such that the omitted terms are analytic near χ = −1, no matter if this point is approached from either of the physical channels χ ∈ (−∞ − i , −1 − i ) or χ ∈ (−1 + i , i ). Thus, the non-analytic behavior of f 1 can be removed, simply by subtracting the antiderivative of the 1 χ+1 · · · term in eq. (6.20), The remainders r 1 and r 2 can only be linear combinations of π 2 log 2 χ, ζ(3) log χ or constants, as they need to be devoid of branch cuts around χ = −1 on all physical sheets (and at any other point other than 0 or infinity). For instance, no dilogarithm nor any term involving log(χ + 1) would have this property. In addition, constraint 2 requires that the second term in eq. (6.16) must be pole-free at χ = −1 and, therefore, we must have that f 2 (−1) = 0 on all physical branches. Imposing this constraint will separately fix the constant term in r 2 and the coefficient of ζ(3) log χ, as the value of log χ is equal to ±iπ, depending on the channel under consideration. It remains to impose the Regge limits. In the case of f 1 , the χ → 0 behavior forces f 1 (0) = 0 hence r 1 = 0. In the case of f 2 , the double-logarithmic term in eq. (6.15) can be used to fix the remaining π 2 log 2 χ ambiguity. The function I 1 (χ) in eq. (6.16) is then uniquely fixed. We find that the function determined this way automatically fulfills the remaining Regge limits given in eq. (6.15), which we view as a nontrivial consistency check.
In conclusion, we have obtained that As a cross-check on this result for I 1 (χ), we have tested it against numerical integration 15 for a number of points with χ > 0 (with ∼ 8 digits precision).
Plugging this into eq. (6.5), we obtain the following complete results for the chiral numerator integrals I ++ and I +− : These formulas are such that, with the standard branch choice for the polylogarithms, the result is real in the Euclidean region χ > 0. Also, we refer to the footnote around eq. (6.1) for an explanation of our conventions. If desired, these results could be rewritten in terms of classical polylogarithms such as Li 4 , but we have not found such rewritings particularly illuminating. Equations (6.23)-(6.24) contain ζ(3) terms which violate the uniform transcendentality degree of the other terms, which may be surprising at first sight. We tentatively attribute this to double-triangle integrals present in the difference between our I +± and the "true" chiral integrals, as discussed in the footnote above eqs. (5.18)- (5.19). It would be interesting to evaluate explicitly the difference and see if the ζ(3)-terms disappear.
As we were hoping, we have found that the chiral numerator integrals I +± admit rather compact analytical expressions. We invite the reader to compare our results against earlier results in the literature for double-box integrals. We refer to eqs. (22)-(25) of ref. [63] for the analytical result for the scalar double box, and to eq. (13) of ref. [64] for the analytical result for the double box with the ( 1 + k 4 ) 2 numerator insertion.

Discussion and conclusions
Generalized unitarity is a method for computing loop-level scattering amplitudes that has been applied very successfully at one loop, in particular to computations of processes with many partons in the final state. In ref. [34], the first steps were taken in developing a fully systematic version of generalized unitarity at two loops. In the approach followed there, the two-loop amplitude is decomposed as a linear combination of basis integrals, in similarity with eq. (1.1). The goal of the calculation is the determination of the integral coefficients as functions of the external momenta; once this is done, the amplitude is determined. At two loops, the integrals with the leading topology in the basis decomposition have the double-box topology, illustrated in figure 1. 16 The integral coefficients of the double-box integrals are determined by applying to both sides of the basis decomposition of the twoloop amplitude so-called augmented heptacuts. These are defined by replacing the seven propagators in the double-box integrand by complex δ-functions. This will freeze seven out of the eight degrees of freedom in the two loop momenta. The Jacobian arising from solving the δ-function constraints contains poles, known as leading singularities, and the remaining integration can be chosen as a contour in the complex plane, enclosing these leading singularities.
Our strategy towards applying the generalized-unitarity method to two-loop QCD amplitudes is not to try to determine the complete QCD integrand, but only the coefficients of a small (possibly minimal) number of craftily chosen "master" integrals, to which all the other ones can be reduced using integral identities (for example, integration-by-parts identities or Gram determinant identities). Projecting out these identities imposes constraints on the contours, namely, that on allowed contours the integral identities which are valid on the physical (uncut) contour must remain valid. As explained in ref. [34], contours satisfying this consistency condition are guaranteed to produce correct results for scattering amplitudes in any gauge theory. 17 These master contours can be chosen so as to extract the coefficient of any particular master integral in a given basis. A perplexing feature of the master contours obtained in ref. [34] is that they are not uniquely defined: indeed, they were found to be characterized by 6 free parameters.
In this paper, we explain this phenomenon as a simple redundancy of variables and find that the two-loop master contours are unique, in perfect analogy with the situation in one-loop generalized unitarity.
Our starting point is a careful examination of the solutions to the heptacut of the double-box integral with an arbitrary number of external states. As mentioned above, the heptacut amounts to setting all the propagators in the double-box graph on-shell. This will freeze all but one of the degrees of freedom in the two loop momenta; for generic external momenta, this degree of freedom is necessarily complex and thus naturally parametrizes a Riemann surface. We have provided a complete classification of the solutions to the heptacut constraints, explained in section 3, based on the number of three-point vertices in the doublebox graph. We find that as the number of three-point vertices is decreased, there are six, four or two Riemann spheres, intersecting pointwise and linked into a chain; or, ultimately, an elliptic curve. We find that the intersection points of these spheres coincide with the poles of the Jacobian arising from linearizing the heptacut constraints; i.e., the two-loop composite leading singularities. In section 5 we explain that this gives rise to identifications which explain the mystery found in ref. [34].
Moreover, we find that at the intersection points of the Riemann spheres, one of the two loop momenta becomes collinear with the massless external momenta attached to the respective vertices of the double-box graph. As integration regions where the loop momentum is becoming collinear with massless external momenta are associated with infrared divergences in the uncut loop integral, this provides a natural physical interpretation of two-loop leading singularities.
In the case when the double-box graph contains no three-point vertices, we find that the Riemann surface associated with the maximal cut is an elliptic curve; i.e., a torus. As discussed in subsection 3.4.1, the presence of an elliptic curve opens up a possibility 17 We stress that the method makes no assumption regarding the powers of loop momentum present in numerators. Thus, as discussed at the beginning of section 5, our method applies indiscriminately to the contributions from planar double boxes in any quantum field theory.
to connect, for the first time precisely, the long-held belief that the maximal cuts of a given integral should be connected to the analytic structure of its integrated expression. In particular, we argued that an integral whose heptacut contains an elliptic curve is unlikely to be expressible in terms of polylogarithms.
The integral coefficients of two-loop amplitudes are functions of the dimensional regulator . The extraction of the O( ) contributions poses a technical difficulty as they may multiply 1 k singularities in the integrated expressions for the two-loop integrals and thus produce finite contributions to the amplitude. Unfortunately, they are not obtainable from cuts performed in strictly four dimensions and must instead be computed by taking cuts D = 4 − 2 dimensions, something which is technically much more involved. One potential way to minimize this technical problem would be to use a basis in which as many integrals as possible are infrared finite. Namely, although the integral coefficients may still have O( ) corrections, these contributions would then multiply infrared finite integrals and hence would not be physically relevant. Of course, as two-loop amplitudes do have infrared divergences, a basis of integrals must necessarily contain some infrared divergent integrals, but it is still plausible that a basis containing as few such integrals as possible will minimize the amount of work needed to obtain their coefficients. In section 5 we have shown that the chiral integrals recently introduced by Arkani-Hamed et al. provide a basis for the double-box contributions to (four-point) two-loop amplitudes in any gauge theory. These integrals are infrared finite, allowing for their coefficients to be obtained from strictly four-dimensional cuts.
In section 6 we evaluated these integrals analytically for the case of four massless external momenta. As hoped, we have found that they are given by very simple and compact expressions (see eqs. (6.23)-(6.24) for our final results). Our computation consisted of a more or less direct Feynman parameter integration, together with a few tricks. It will be very interesting to see if a similar finite basis of masters is available in the case of five particles, in the double-box and pentagon-box topologies. Moreover, we are very hopeful that such finite integrals can be computed analytically by some means.
Using the integral identities which have been explicitly constructed at five points [54], we believe it will be possible to explicitly construct the (unique) master contours in this case using our results. We hope to return to this question in the near future. In general, a better understanding of integration-by-parts and Gram determinant identities would be very helpful. Although in this work we have not considered other topologies with seven or more propagators, such as pentagon-boxes, we expect the extraction of their coefficients to be much simpler than that of the double boxes, due to the presence of explicit octacuts. It will also be important to understand their extraction in the future. In addition, it is important to understand the extraction of coefficients of integrals with a smaller number of propagators. Other natural extensions of our work would include non-planar topologies, or integrals with internal masses. We are hopeful that the ideas presented in this paper will be useful for answering these questions.
Another important question concerns the dimension of the space of master integrals of a given topology. Although for this particular question we have only considered the case of four particles, we conjecture that in the general case, the master contours are still unique, and that their number is precisely equal to the number of independent master integrals of the corresponding double-box topology. Verifying this would be very interesting.

A Residues of maximally cut amplitudes
An amusing application of the identities (3.15)- (3.20) arises in the context of taking residues at Jacobian poles of the heptacut two-loop amplitude J(z) 6 j=1 A tree j (z) S i . In generalizedunitarity calculations, this quantity forms the input out of which the integral coefficients of the two-loop amplitude are computed. The identities (3.15)-(3.20) relate, in particular, the residues of the heptacut two-loop amplitude across different kinematical solutions. As a result, they explain the vanishing of certain residues, as well as the seemingly accidental equality between pairs of other residues. This in turn allows one to cut the work of evaluating these residues in half.
As an example, let us consider the heptacut illustrated in figure 9 of the two-loop amplitude A (2) (1 − , 2 − , 3 + , 4 + , 5 + ). The helicities assigned to the internal and external states allow only gluons to propagate in the loops; this in turn implies [65] that the results for the heptacut amplitude within N = 4, 2, 1, 0 Yang-Mills theory are identical.
The heptacut J(z) 6 j=1 A tree j (z) S i shown in figure 9 of this amplitude only receives nonvanishing contributions from kinematical solutions consistent with the internal helicites shown there -in particular kinematical solutions whose (3, p 1 , p 2 )-vertex is MHV. By inspection of figure 3 we thus see that For the heptacut amplitude evaluated at the remaining three kinematical solutions, direct One observes that the residues at z = 0 in solutions S 4 and S 5 are vanishing. This can now be easily explained by eqs. (3.19) and (3.16) as a consequence of the vanishing (A.1) of the heptacut amplitude on solutions S 6 and S 2 , respectively. Moreover, eq. (3.18) relates the residues at the nonzero Jacobian poles in S 3 and S 4 ; similarly, eq. (3.17) relates residues in S 3 and S 5 . In conclusion, the identities (3.15)-(3.20) allow us to cut the work of evaluating the residues of the heptacut two-loop amplitude at the Jacobian poles in half.

B An elliptic curve in planar N = 4 super Yang-Mills
Computations in N = 4 super Yang-Mills (SYM) theory are generally much easier than in QCD or pure Yang-Mills, in no small part due to its so-called dual conformal (super)symmetry, which is a hidden symmetry of the planar limit of the theory, invisible from its Lagrangian [66]. Several two-loop amplitudes have been computed analytically so far, and all have been expressible in terms of special functions called degree-4 (multiple) polylogarithms. In some cases, using this information, it was even possible to guess nontrivial amplitudes [67,68]. Therefore, it is an important question whether all amplitudes in planar N = 4 SYM are given by polylogarithms.
Here we wish to give evidence that N = 4 SYM knows about much more than polylogarithms, even when only massless internal states are present. We will do so by exhibiting a specific helicity configuration for 10 particles which will be given by the single integral in figure 6, plus nothing else, ruling out any possible cancelations. As argued in the main text, we find it extremely unlikely that this integral is expressible in terms of polylogarithms.
In ref. [74], the super Wilson loop was expressed in the form are vertex factors, and the edges are E i = Pe −i Aµdx µ +O(χ) . The k i are the on-shell momenta of the scattering amplitudes, which are also the lengths of the Wilson polygon segments. As we will now argue, the chosen Grassmannian component is cooked up to make the edge factors 1 at this order; that is, E i → 1. To see this, we will show that for this component there is no coupling to the fermions of the theory. This will rely on very general properties of the edge interactions, which is why we do not give explicit expressions here. For instance, the edge E i contains a coupling of the form ψχ i χ i χ i , which vanishes here because no cubic power of any given χ i is present in the Grassmann component under consideration. The vertex factor between k i and k i+1 contains cubic interactions of the form χ i χ i χ i+1 , but no such product can form a 4 of SU(4) R for the chosen component. We conclude that all the couplings of the Wilson loop to ψ vanish. Direct couplings to the field strength tensor F µν vanish for the same reasons. Regarding the couplings toψχ, these come from edge integrals of the form E ∝ ψ χ. If such a coupling were to contribute, theψ propagator would have to land on a ψ field, which could only come from a Lagrangian interaction term of the form g YM ψψφ. While this interaction exists, it is too inefficient a way to use one power of g YM , for an amplitude which must couple 12 essentially distinct χ's using only four powers of g YM (it could only contribute at order g 6 YM ). We conclude that with the chosen component, the Wilson loop couples only to the scalars of the theory. On the edge corresponding to k i there is a φχ i χ i coupling, which could in principle contribute for edges 5 and 10 because our component does contain two χ 5 and two χ 10 . But this would require the next two χ's, for instance χ 4 6 and χ 4 7 , to also couple to a scalar. However, this would not be allowed by the SU(4) R symmetry. We conclude that the edge factors E i are trivial for the chosen Grassmannian component, and that only the vertex factors shown in eq. (B.2) contribute.
Thus, the expectation value of the supersymmetric Wilson loop, for the chosen component, reduces to a correlation function of 6 scalars, each of which couples to two χ Grassmann components. The six scalars in turn propagate to two g 2 YM [φ, φ] 2 interactions from the Lagrangian, producing a contribution of order g 4 YM . All possible distributions of the scalars among the vertices have to be considered, but for the chosen component it is possible to see that only one distribution is allowed by the SU(4) R symmetry. It gives rise to the Feynman diagram shown in figure 10; this is the one and only diagram which contributes. The seven scalar propagators in the diagram reproduce precisely the seven denominators in the right hand side of eq. (B.1) -the Feynman graphs are dual to each other. That is, the momentum space expression for one is equal to the configuration space expression for the other. We omit the verification that the prefactor works out as claimed.
This agreement concerns the integrand in exactly four dimensions. However, the chosen component is infrared finite at two loops because, by infrared exponentiation theorems, any IR divergence of the two-loop amplitude would have to be canceled by IR divergences in the tree and one-loop amplitudes. But as the latter vanish (in D dimensions) for this process, the two-loop amplitude is necessarily IR finite. For this reason, it is certainly enough to consider only the four-dimensional integrand in this case. We have also checked the present result for the four-dimensional integrand against a computer implementation of the recursion relation for the integrand [75], finding perfect agreement.