Differential equations on unitarity cut surfaces

We reformulate differential equations (DEs) for Feynman integrals to avoid doubled propagators in intermediate steps. External momentum derivatives are dressed with loop momentum derivatives to form tangent vectors to unitarity cut surfaces, in a way inspired by unitarity-compatible IBP reduction. For the one-loop box, our method directly produces the final DEs without any integration-by-parts reduction. We further illustrate the method by deriving maximal-cut level differential equations for two-loop nonplanar five-point integrals, whose exact expressions are yet unknown. We speed up the computation using finite field techniques and rational function reconstruction.

After IBP reduction to a basis of master integrals, the master integrals still need to be evaluated, e.g. using the method of differential equations (DEs) [25][26][27][28][29][30]. A recent breakthrough was Henn's canonical form of DEs [31,32], allowing a large class of loop integrals to be expressed in terms of iterated integrals of uniform transcendentality. Various algorithms and software packages [32][33][34][35][36][37][38] have appeared to find algebraic transformations of DEs to the canonical form, while a complementary approach is finding master integrals in the d log form with unit leading singularities [32,39,40].
In constructing DEs, an important intermediate step is IBP reduction, which brings the RHS of the DEs into a linear combination of master integrals. External momentum derivatives increase the power of propagator denominators, so unitarity-compatible IBP reduction is not directly applicable. To solve this problem, one approach is to decrease the power of propagator denominators using dimension shifting [22]. However, we propose an alternative approach that completely avoids doubled propagators, even in intermediate steps. We promote unitarity cut surfaces to be objects embedded in the space of not only loop momenta, but also external momenta. By combining external and loop momentum derivatives to form tangent vectors to unitarity cut surfaces, doubled propagators cancel out, in direct analogy with unitarity-compatible IBP reduction.
It has been proposed that the maximal cut can provide valuable information about differential equations and the function space of the integrals [41][42][43]. The latter two references rely on consistent definitions of unitarity cuts in the presence doubled propagators (see also [44,45]). This issue is bypassed in our approach by construction. Section 2 sketches the basics of our formalism. Section 3 uses inverse propagator coordinates, also known as the Baikov representation in the d-dimensional case, to present the detailed formalism through the one-loop box example. Section 4 applies the formalism to the nonplanar pentabox at the maximal cut level, and finds a system proportional to the dimensional regularization parameter , for tensor integrals with unit leading singularities. Section 5 discusses the use of finite field techniques and rational function construction in speeding up the nonplanar pentabox computation. Some concluding remarks are given in Section 6.
We will explore a similar unitarity-compatible approach to differential equations, with no Feynman integrals involving doubled propagators appearing on the LHS or RHS of the differential equations, even before IBP reduction is performed to simplify the RHS. The advantage is two-fold. First, the differential equations may be put on unitarity cuts, and therefore can be constructed by merging incomplete results on a spanning set of unitarity cuts. Second, unitarity-compatible IBP reduction can be used to reduce the RHS into the original set of master integrals, since no integrals with doubled propagators are ever generated.
We wish to compute the derivative of a Feynman integral with propagators 1/z j and a tensor numerator N , where β µ i is a linear combination of external momenta p j , We are free to add total divergences to Eq. (2.1) without changing its value after integration, obtaining In the above expressions, v µ = v µ (l, p) has polynomial dependence on internal and external momenta, with one free Lorentz index. The final expression Eq. (2.4) has no doubled propagators ∼ 1/z 2 j , if the following condition is satisfied, where f j has polynomial dependence on internal and external momenta.

Relation to unitarity cut surfaces
As we will see, the condition Eq. (2.5) has a nice geometric interpretation in terms of unitarity cut surfaces, similar to what was observed [20] in unitarity-compatible IBP reduction. Consider a L-loop Feynman diagram topology with N internal propagators. For a subset ∆ of all inverse propagators 1, 2, . . . , N , the unitarity cut surface labeled by ∆ is defined as the hypersurface of all points (l 1 , l 2 , . . . , l L ) in the complex loop momentum space which solves the generalized unitarity cut condition, (2.6) Notice that generalized unitarity cuts differ from traditional unitarity cuts in QFT textbooks, treated by e.g. the optical theorem, Cutkosky rules [46], and the largest time equation [47,48], in several respects, 1. Complex rather than real loop momentum space is considered. In particular, the energy component of the cut loop momentum is not required to be (real) positive. The algebraic closure of the complex field guarantees that the unitarity cut condition has solutions for generic external momenta, as long as not too many propagators are cut (which will be assumed to be the case for the rest of the paper). 1 2. There is no direct connection with the discontinuity of Feynman diagrams across branch cuts.
When ∆ contains all inverse propagators, i.e. when the unitarity cut condition sets all the propagators on-shell, the unitarity cut is called a maximal cut. Many problems, such as integrand construction and integration by parts, are simplest at the level of the maximal cut. When a proper subset of propagators are set on-shell, we have a non-maximal cut.
In general, for a diagram topology with N propagators, there are 2 N different unitarity cut surfaces, since we may choose each propagator to be either cut or not cut. Since this paper is concerned with differential equations w.r.t. external momenta, we define "extended unitarity cut surfaces" which is embedded in the space of not only loop momenta but also external momenta. Here the loop momenta are still allowed to be complex but external momenta are required to be real. As usual, the surface is defined by the cut condition, Eq. (2.6).
For any of the 2 N extended unitarity cut surfaces ∆ with z k ∈ ∆, the cut condition sets z k = 0, while the condition for the absence of doubled propagators, Eq. (2.5), becomes Therefore, the expression which we refer to as a "DE vector", is a tangent vector to every extended unitarity cut surface embedded in the space of both internal and external momenta. The loop part, v µ ∂ ∂l µ , (2.9) is called an "IBP vector". In the special case that the DE vector has no external momentum derivative, the IBP vector itself is a tangent vector to unitarity cut surfaces, and is used in unitarity-compatible IBP reduction. Computational algebraic geometry is used to find IBP vectors in the literature, and will also be used in this paper to find DE vectors.
3 Detailed formalism in inverse propagator coordinates

Inverse propagator coordinates
We re-examine the simple example of the one-loop box [26] using our method. We assume that all internal and external lines are massless. The scalar box integral with some numerator N , shown in the leftmost diagram of Fig. 1, is with The kinematic invariants are It is well known that after IBP reduction, there are 3 master integrals. For our purposes, they are conveniently chosen as the scalar box I box , the s-channel triangle I The loop integral can be parameterized in the inverse propagator coordinates, in either 4 or d dimensions. This parameterization goes back to Cutkosky's proof of the cutting rules [46], and has been systematically studied by Baikov for the d-dimensional case [49][50][51]. More recently, this parameterization was used in [20,21] for unitarity-compatible IBP reduction. A detailed explanation of the Baikov representation recently appeared in [43], where a public code was made available, and applications to differential equations were discussed. We now derive this parameterization for the one-loop box, and refer the readers to the literature for the multi-loop case.
Using the Van Neerven-Vermaseren basis [52], the metric tensor η µν is written as the sum of a "physical" component in the 3-dimensional space spanned by external momenta, and a "transverse" componentη µν in the remaining (d − 3)-dimensional space which is orthogonal to every external momentum, where we used the inverse of the Gram matrix G, We will need the relations, In the above equations we have defined where F is called the Baikov polynomial. For the one-loop box, it evaluates to The loop integration measure, multiplied by the propagators, is re-written as In the second-to-last line of the above equations, we integrated outμ 2 against the delta function, and in the last line we made a linear transformation of integration variables with unit Jacobian. Eq. (3.19) accomplishes the transformation of the loop momentum from the Lorentzian / Euclidean component coordinates to the inverse propagator coordinates. We will adopt the differential form notation for IBP relations [20,21], which helps to transform between different coordinate systems for the loop momenta. We write where Ω is a maximal differential form in the loop momentum space. For convenience, we define Ω to include both the integration measure and the propagators. The total divergence integral of any

Vectors for differential equations
We first discuss IBP vectors, i.e. the loop component of DE vectors. The IBP vector Eq. (2.9) can be transformed into the new parameterization, Since v µ is required to have polynomial dependence on external and internal momenta, it is easy to show that in the new coordinates, v i must have polynomial dependence on the inverse propagators z. Furthermore, v µ ∂ µ is required to satisfy rotational invariance in the where the proportionality constant has polynomial dependence on z i . In the inverse propagator coordinates, this means for some v F which is a polynomial in the z variable. 3 The above equation is a crucial criterion for a valid IBP vector, derived in a different way in [21].
with v F defined in Eq. (3.25). Setting the tensor numerator N to 1 in Eq. (2.3), the 2nd term in the square bracket is given by Eq. (3.26), while the first term is simply, suppressing the i index, Adding (3.27) and the vanishing integral Eq. (3.26), we obtain, The above expression has no doubled propagators when for some polynomial f k = f k (z) for every inverse propagator z k . This is equivalent to Eq. (2.7). In summary, in inverse propagator coordinates, a unitarity-compatible DE vector that does not lead to doubled propagators is required to satisfy two requirements, Eq. (3.25) and (3.31). Combining these two equations gives Given the derivative against external momenta, β µ ∂ µ , Eq. (3.33) is a polynomial equation in the unknown polynomials f j and v F . This is an inhomogeneous version of the "syzygy equation" of Ref. [21]. 4 We use the computational algebraic geometry package SINGULAR [53] to find a particular solution for f j and v F in Eq. (3.33), which in turn fixes v k in Eq. (3.31). This allows us to compute the derivative of the loop integral w.r.t. external momenta using Eq. (3.29).

Evaluating DEs: one-loop box warm-up
Consider the derivative with respect to t = (p 2 + p 3 ) 2 , which annihilates s = (p 1 + p 2 ) 2 and keeps external momenta on-shell. Combined with Eqs. (3.2)-(3.5), we find the following expressions for the RHS of Eq. (3.33), Solving Eq. (3.33) for the unknown polynomials f i and v F , we find the following solution, which can be easily checked, where we set d = 4 − 2 . The triangle integrals are one-scale integrals whose derivatives against t can be fixed by simple dimensional analysis, so we omit the calculation. Transforming to the basis suggested in Refs. [31,32], we obtain which agrees with the result obtained from the standard approach involving IBP reduction of integrals with doubled propagators. In this simple example, we directly obtain the DEs without any IBP reduction, but unitarity-compatible IBP reduction of tensor integrals is needed in more complicated cases.

4D leading singularities
Although this paper mainly concerns unitarity cuts in d dimensions, the inverse propagator coordinates also make it easy to compute cuts in 4 dimensions. Setting d = 4 in Eq. (3.19), the maximal cut residue of the box integral I box is equal to where F (0) = F (z i = 0). Therefore, st I box is an integral whose leading singularity is a Q number with no dependence on s or t. This is an important criterion for selecting candidate integrals with uniform transcendentality. We now translate the pentagon integrals with unit leading singularities [39] to the inverse propagator coordinates, which serves as a building block for the master integral choice for the nonplanar pentabox. This loop topology is shown in Fig. 2. The five kinematic invariants can be chosen, in a cyclic invariant fashion, to be s 12 , s 23 , s 34 , s 45 , s 51 , where s ij = (p i + p j ) 2 . The Gram matrix is again defined as G ij = p i · p j , 1 ≤ i, j ≤ 4. In D dimensions, the pentagon integral with a numerator N can be expressed in terms of the five inverse propagator coordinates, z i = (l − i−1 j=1 p j ) 2 , whereas in 4 dimensions, the Baikov polynomial resides in a Dirac delta function since the z i 's become linearly dependent [20],  Choosing then Eq. (3.43) produces the leading singularity 1 on the cut solution z 5 = z − 5 , and 0 on the other cut solution z 5 = z + 5 . One important feature which makes Eq. (3.46) a valid ansatz is that the z 5 -independent term, −F (0)/ √ det G, is a cyclic invariant expression which is the same on all the five possible 4D maximal cuts. So Eq. (3.46) can be promoted to the full uncut expression for N by combining the z-dependent terms obtained on the individual cuts, and the result corresponds to the chiral pentagon integral of Ref. [39] up to terms that vanish on 4-dimensional box cuts.

A non-planar five-point topology at two-loops
To illustrate our method for computing differential equations on unitarity cuts, we consider two-loop five-point scattering. While planar master integrals have been computed [54,55], the nonplanar counterpart remains unknown. Ref. [43] suggested probing the properties of such integrals via the maximal cut. We focus on the nonplanar pentabox in Fig. 3, and compute differential equations at the maximal cut level. The inverse propagators are while the irreducible numerators are initially chosen as Using inverse propagator coordinates, the loop integral can be written as, omitting overall factors, where det G is the Gram determinant for the external momenta, defined in the same way as for the one-loop pentagon, and F (z) is defined by the (d − 4)-dimensional components of l 1 and l 2 , which are µ 1 and µ 2 , by The five independent kinematic invariants are defined in the same way as for the one-loop pentagon, as s ij = (p i + p j ) 2 , (ij) = (12), (23), (34), (45), (51) . We perform IBP reduction at the maximal level for tensor numerators with up to 6 total powers of z 9 , z 10 and z 11 , using an in-house implementation of the algorithm of [21]. This involves solving a syzygy equation, which is the same as our Eq. (3.33) with the RHS set to zero, to find vectors whose total divergences generate IBP relations without doubled propagators. Three master integrals, I 1 , I 2 , and I 3 , are found, with tensor numerators respectively. Next, the derivatives w.r.t. the four scaleless kinematic invariants, ∂/∂χ ij , are each written in terms of external momentum derivatives ∼ β µ (∂/∂p µ ), in the fashion of Eq. (3.34). This fixes the RHS of Eq. (3.33). On a unitarity cut, Eq. (3.33) is simplified, allowing a solution to be found by the computer more quickly. 5 Finding one particular solution for f j and v F using SINGULAR gives us the needed ingredients to compute the derivative of the loop integral in Eq. (2.4). We obtain a combination of tensor integrals without doubled propagators, which are then reduced to the three master integrals using the IBP reduction procedure described above. The end results are the differential equations relating the three master integrals with their partial derivatives against the four scaleless kinematic invariants, The full results for the matrix M ij AB are attached in an ancillary file pentaCrossBox.m in the Mathematica format. We will elaborate on the computation techniques in Subsection 5.2. A sample component is which have poles that can be identified with kinematic singularities, such as χ 45 ∝ (p 4 +p 5 ) 2 , and (χ 23 − χ 45 + 1) ∝ (p 1 + p 3 ) 2 . Since Eq. (4.8) allows second derivatives to be computed, we are able to perform another check using the consistency condition [35], As is the case for the component shown in Eq. (4.9), the complete matrix is linear in but not proportional to , because we have not yet transformed the result into a "good" basis of master integrals with unit leading singularities. Such a "good" basis can be found easily using the 4D cut method in [31]. Roughly speaking, cutting the l 2 box sub-loop produces the Jacobian 1 z 5 z 5 , where z 5 ≡ (l 1 + p 5 ) 2 , z 5 ≡ (l 1 + p 4 ) 2 . (4.11) So we can recycle the one-loop chiral pentagon expression Eq. (3.46) to write down three tensor integrals,Ĩ 1 ,Ĩ 2 , andĨ 3 with unit leading singularities. Their numerators are (ignoring the dependence of the chiral pentagon numerator N chiral on z 1 , z 2 , z 3 , z 4 which vanish on the maximal cut of the nonplanar pentabox), respectively. The first two of these integrals are among the nonplanar N = 4 SYM integrands for this particular topology given in [40].
The new basis is related to the old one via a matrix T , 13) and the differential equations in the new basis, (4.14) are related to the old system, Eq. (4.8), by the transformation formula, After the transformation, the matrices are proportional to = −(d − 4)/2, and now involve not only polynomials but also square roots of the Gram determinant. The square roots are eliminated by switching to momentum-twistor variables [56] using the parameterization given in Appendix (A.2) of Ref. [57]. The momentum-twistor variables, x i with i = 1, 2, 3, 4, 5, are related to the usual kinematic invariants by where tr 5 is defined via the Gram determinant, Since there are only 4 dimensionless ratios of kinematic invariants, we will fix x 1 = s 12 = 1, effectively only looking at the dependence of the integrals on x 2 , x 3 , x 4 , x 5 . After re-writing the differential equations in terms of the above momentum-twistor variables, we used the CANONICA software package [36] to transform the differential equation into a form that contains "dlogs" [31] (a few seconds of computation time is used), In the above equation, I is a column vector consisting of the 3 maximal-cut master integrals of unit leading singularities,Ĩ 1 ,Ĩ 2 , andĨ 3 , defined via the tensor numerators in Eq.  19) and the explicit expressions for the 3 × 3 matrices are, This form of the maximal-cut differential equations indicate that the solutions are multiple polylogarithms involving the 11 symbol letters in Eq. (4.19), with uniform transcendentality in the expansion, as explained in Refs. [31,32].

Double box and comparison with FIRE5
Having illustrated our method at the one-loop level in Section 3 and at the two-loop level in Section 4, we test the method for the double box and compare with existing methods. Being not too simple or too complicated, the double box topology can be handled by both traditional and unitarity-based methods, allowing for a meaningful comparison. The double box with massless internal and external lines is shown in Fig. 4, with kinematic variables defined in the caption. It is well known [58] that there are 2 top-level master integrals, which may be chosen as the scalar integral I s dbox and the tensor integral With daughter topologies included, there are 12 master integrals in total. Among them, 8 master integrals are independent after accounting for the discrete symmetry, which leaves the 7 propagators invariant up to a permutation. We re-compute the most non-trivial part of the system of differential equations, i.e. the t-derivative of the 2 top-level master integrals, expressed in terms of the 8 master integrals. Our method is applied to a spanning set of 6 different unitarity cuts, The discrete symmetry Eq. (5.2) is used to speed up the calculation by relating the cut Γ 2a with Γ 2b , and Γ 3a with Γ 3b . In particular, the IBP relations generated on one cut automatically become valid IBP relations on the related cut after the symmetry transformation. Not surprisingly, these are essentially the same cuts used in Ref. [21] for IBP reduction of double box integrals. By merging the results on the 6 cuts, we reproduce the following differential equations, where daughter topology integrals are fully computed but omitted in the above sample results. Both FIRE5 and our own unitarity-based code begin with "preparation runs" purely in Mathematica to process the diagram topology information supplied by the user. For the "final runs", FIRE5 is used in the C++ mode, while our own code is written in Mathematica and SINGULAR. Only the final runs, which reflect the true computational complexities, are included in our timing comparison. The time required to obtain the results in Eq. (5.9) is shown in Table 1.
Though there could be more improvements by switching to statically compiled computer languages such as C++, our code already offers significant improvement in speed over a calculation based on FIRE5, for the following possible reasons, 1. The lack of doubled propagators reduces the number of different integrals that appear in IBP relations, because for combinatorial reasons, the number of integrals grow rapidly when propagators are allowed to be raised to higher powers.
2. The use of unitarity cuts reduces the computational complexity.
3. Besides solving syzygy equations, SINGULAR is also used to solve the sparse linear system formed by the IBP relations. Unfortunately, due to lack of a controlled comparison, we do not know how this affects performance compared to Fermat [59] used in the C++ version of FIRE5.

Finite field techniques and rational function reconstruction
The computation of the differential equations for the nonplanar pentabox posed challenges in terms of CPU time and memory consumption. Given that the results, e.g. Eq. (4.9), are rational functions in χ ij , we use the rational function reconstruction technique of [60] to fit the analytic result from numerical inputs of χ ij . The algorithm of [60] reconstructs multivariate rational functions in two steps, (i) fitting univariate rational functions, and (ii) fitting multivariate polynomials, using the input from many iterations of step (i). We use a simple private implementation which performs most of the work by exploiting the built-in capabilities of Wolfram Mathematica. 6 6 In Mathematica 10, step (i) is accomplished by the command FindSequenceFunction with the option FunctionSpace -> "RationalFunction".
Step (ii) is accomplished by the command InterpolatingPolynomial. The latter command allows the option of computing in a finite field Zp.
For the nonplanar pentabox computation, step (i) requires 18 kinematic points for each iteration, and 495 iterations are performed to produce the input for step (ii).
Step (ii) fits polynomials in 3 variables, with degrees up to 8. Finite field techniques are used to accelerate step (ii): using a large prime p, the full result can be constructed from its image in Z p probabilistically using a minor modification of the extended Euclid algorithm [60,61]. After completing steps (i) and (ii), the fitted results are validated against new computations with additional random rational values of χ ij .
Here we give more information to quantify the performance gains from rational function reconstruction and finite field techniques. The computation of the differential equations for the nonplanar pentabox, with analytic dependence on the kinematic invariants, is very time consuming and does not finish after 48 hours. 7 However, with (rational) numerical kinematic invariants, the computation finishes in a few seconds per kinematic point on a modern computer. A total of 28 hours is used in evaluating differential equations on 8910 kinematic points and reconstructing the full analytic results. The last step of the calculation, i.e. multivariate polynomial fitting, dramatically benefits from finite field techniques and takes about 52 seconds to reconstruct all the 36 entries of the four 3 × 3 matrices. In contrast, when finite field techniques are turned off, about 395 seconds are needed to reconstruct only one of these 36 entries (with the computation aborted afterwards), which is slower by more than 2 orders of magnitude.

Conclusions
We have proposed a new method for constructing differential equations for Feynman integrals, which avoids generating integrals with doubled propagators, instead producing tensor integrals to be reduced by unitarity-compatible IBP reduction. In fact, for the simplest cases such as the one-loop box, no IBP reduction is needed at all. Our method allows constructing differential equations from a spanning set of unitarity cuts in d dimensions, with IBP reduction also performed on the cuts.
Applying our method to the nonplanar pentabox, we obtained the homogeneous differential equations on the maximal cut in Henn's canonical form. This allows us to confirm that the master integrals, at least when evaluated on the maximal cut, are multiple polylogarithms with uniform transcendentality in the expansion. We have extracted the 11 symbol letters, which are polynomials of momentum-twistor variables.
We also demonstrated that finite field techniques and rational function reconstruction, which are emerging as new tools in studying scattering amplitudes [60][61][62], are useful in computing differential equations for Feynman integrals.
There are several possible directions for follow-up studies. One direction is extending the calculation to other nonplanar five-point topologies, which can be done straightforwardly. It would be desirable to construct an automated implementation of our method, perhaps as an extension to unitarity-compatible IBP reduction software packages, such as Azurite [22], since many computation steps can be shared. Eventually, we would like to construct full DEs for nonplanar two-loop five-point integrals, which are relevant for NNLO QCD corrections for 2 → 3 scattering processes at the LHC [54,55].