Maximal Cuts in Arbitrary Dimension

We develop a systematic procedure for computing maximal unitarity cuts of multiloop Feynman integrals in arbitrary dimension. Our approach is based on the Baikov representation in which the structure of the cuts is particularly simple. We examine several planar and nonplanar integral topologies and demonstrate that the maximal cut inherits IBPs and dimension shift identities satisfied by the uncut integral. Furthermore, for the examples we calculated, we find that the maximal cut functions from different allowed regions, form the Wronskian matrix of the differential equations on the maximal cut.


Introduction
Quantum field theory scattering amplitudes are mathematical quantities enabling physicists to make predictions for physical observables in high energy particle experiments such as the Large Hadron Collider (LHC) at CERN. Although scattering amplitudes are among the most important objects in this research direction, sufficient precision frequently requires explicit computations that are extremely challenging, even with powerful modern techniques. The reason is often the complexity of the Feynman integrals involved, and an inadequate understanding of the underlying mathematics and the surprisingly rich hidden structures of scattering amplitudes that are continuously being unravelled.
The past two decades have seen enormous progress in the development of new enhanced methods for computing multiloop scattering amplitudes. The traditional techniques due to Feynman are no longer preferred by experts for state-of-the-art calculations. Instead, scattering amplitudes are typically reduced to a linear combination of integrand or integral basis elements, whose coefficients then become the primary quantities of interest after the integrals have been carried out once and for all. All one-loop integrals can be expressed in terms of simple algebraic functions along with the logarithm and dilogarithm, whose arguments are again algebraic functions. Today, fully automated computation of one-loop amplitudes has been achieved, either via the unitarity method [1,2] and its refinements [3,4], or by the Ossola, Papadopoulos and Pittau (OPP) approach [5,6] at the level of the integrand. More recently, extensions of these techniques to two loops in general theories have been reported, forming the frontier of next-to-next-to-leading (NNLO) corrections. A key element in these developments has been the application of (computational) algebraic geometry [7]. See refs. [8][9][10][11][12][13][14][15] for the multiloop version of the OPP method, and refs. [16][17][18][19][20][21][22][23][24][25][26] for progress on direct extraction of integral coefficients from an integral basis.
One of the remaining bottlenecks in the unitarity and OPP based methods at the multiloop level is the computation of the Feynmans integrals themselves. Morover, at two loops and beyond it is more complicated to even determine an appropriate integral basis [27] at higher multiplicity. Given a complete set of master integrals for the problem in consideration, the standard procedure for evaluating them is to derive differential equations [28][29][30][31][32][33] in the external kinematic invariants, reduce the resulting expression using integration-by-parts (IBP) identities [34] to form a linear system of equations. The integrated expressions are constructed from a much less restricted class of transcendental functions, including for instance generalized polylogarithms. These ideas have proven extremely useful in practice over the years. In particular, if the basis integrals are chosen properly, the differential equations are brought to the canonical form ( -form) proposed by Henn [35], leading to significant simplifcations. (See ref. [36] for a pedagogical review and refs. [37][38][39][40][41][42] for algorithms and packages for finding the canonical form.) Motivated by the tremendous success of generalized unitarity at one loop, we present here a systematic strategy for evaluating maximal cuts of multiloop Feynman integrals, properly defined within dimensional regularization. D = 4 two-loop maximal unitarity was first achieved by the elegant contour method in ref. [16] by Kosower and Larsen, and then generalized to other 4D integrals with external or internal massive legs, nonplanar topology and three-loop order [17][18][19][20][21][22][23][24]. The work by Frellesvig and Papadopoulos [43] studied the D-dimensional maximal cut via the Baikov representation [44][45][46], and also explicitly presented the expansion around 4D for the maximal cut function. On the other hand, using the multivariate residue calculus of Leray, refs. [47,48] give a precise definition of the cut Feynman integral in dimensional regularization, and argue that integral relations carry over from the uncut Feynman integrals to the cut integrals. While refs. [47,48] focus on the one-loop case, they predict the construction works at higher loops as well.
In this paper, we systematically study maximal cuts in any spacetime dimensions, by computing the Baikov integrals on the maximal cut over all possible regions, and verify that the maximal cut functions automatically incorporate all integral relations such as IBPs, dimension recurrence relations, and differential equations on the maximal cut. This method applies equally well to planar and nonplanar integrals with and without massive particles. A careful treatment lends credence to the belief that for an integral topology with m master integrals, each master integral would have precisely m linearly independent maximal cut functions in D dimensions. We provide nontrivial evidence that these m maximal cut functions form the Wronskian matrix associated with the differential equation satisfied by the master integrals on the maximal cut. The leading terms of this Wronskian matrix are useful to transform the differential equation to the canonical form. This paper is related to ref. [49] by Primo and Tancredi, and ref. [50] by Zeng, which study the differential equations on the maximal cut systematically, and ref. [43] by Frellesvig and Papadopoulos which applies the efficient loop-induction method to find the maximal cut in the Baikov representation. We remark that our paper is characterized by (i) always retaining complete dimension dependence for all cuts in closed form, so the limit behaviour of the cut function near any integer dimension can be easily obtained, (ii) giving full dependence of the irreducible scalar product (ISP) indices to make all integral relations (like IBP relations) manifest, (iii) most importantly, providing the complete solution system (Wronskian matrix) for D-dimensional IBPs, dimension recurrence identities and differential equations on the maximal cut, from an analysis of all allowed integration regions of the Baikov representation on the cut. This paper is organized as follows: in section 2, we present the Baikov integral representation with the maximal cut and integral regions for real kinematics. In section 3 and 4, we review the simple D-dimensional maximal cut examples with zero or one ISP. Section 5 and 6 contain our main examples, for which the Baikov integration on the cut over different regions gives independent solutions for on-shell IBPs, dimension recurrence relations and differential equations. We explicitly show these solutions are complete by studying the Wronskian of the differential equation.

Baikov representation and maximal cuts
We are interested in L-loop Feynman integrals with n external momenta and k propagators, where m = (n − 1)L + L(L + 1)/2 and k ≤ m. D 1 , . . . , D k are denominators of Feynman propagators, and D k+1 , . . . , D m are the irreducible scalar products (ISPs). So we require that (for integrals in this particular sector), We use the Baikov representation [44][45][46] of (2.1). Schemetically, where F (z) is the Baikov polynomial. The overall factor C(D, x) is a product of hypersphere areas, the Jacobian of the Baikov transformation and the Gram determinant. The kinematic variables are collectively called x.
In this paper, we simply consider real-valued external and internal momenta to simplify the discussion of the Baikov integration region A. For real momenta, A is determined by the spacetime metric signature and Cauchy-Schwarz inequality. For L = 1, the integration region A is simply defined by F (z) ≥ 0. For L = 2, the integration region A is defined by F (z) ≥ 0, µ 11 (z) ≥ 0 and µ 22 (z) ≥ 0, where µ 11 (z) and µ 22 (z) are defined as following: the loop momenta are separated into the projections in the (n − 1)-dimensional space (spanned by external momenta) and the orthogonal complement, The inner products of l ⊥ 1 and l ⊥ 2 are . In terms of the Baikov representation, the µ's become polynomials in z's. For real internal momenta, µ 11 ≥ 0, µ 22 ≥ 0. Furthermore, by the Cauchy-Schwarz inequality, Unitarity cuts become manifest in the Baikov representation. For example, the maximal cut in Baikov representation is to consider the multivariate residue at z 1 = z 2 = . . . = z k = 0 [43][44][45][46][51][52][53]. If the integral has no double propagators, i.e., a 1 = a 2 = . . . = a k = 1, the maximal cut becomes where Ω is the intersection of A and the hyperplane z 1 = z 2 = . . . = z k = 0. For the case with some a i > 1, 1 ≤ i ≤ k, derivates of the Baikov polynomial are needed to get the residue. This form can be used to derive integration-by-parts (IBP) identities [52,53] and differential equations [50] on the maximal cut, and to identify master integrals [54,55], by using Morse theory, tangent vectors and syzygy computations. More generally, the non-maximal cut in Baikov form can be used to derive the complete set of IBPs [52,53].
In this paper, we systematically study (2.5) in detail. We find that frequently the region Ω decomposes into several subregions, (see figure 5 for the subregions of massless double box as an explicit example) where on the boundary ∂Ω j of each subregion Ω j , F = 0. We denote by s the number of such subregions. Then we can explicitly carry out the integration on each Ω j and apply analytic continuation in x and D. The resulting function is named as the maximal cut function on the subregion Ω j , Since F = 0 on ∂Ω j , the possible surface term from the integration of a total derivative vanishes. Hence it is clear that for each fixed j, the functions I 1,...1,a k+1 ,...an (j) m.c. satisfy (the same form of) integration-by-parts identities on the maximal cut. Similarly, for each fixed j, I 1,...1,a k+1 ,...an (j) m.c. 's satisfy (the same form of) dimension shift identities and differential equations on the maximal cut.
The integrals over these subregions may not be independent. For each j, we may consider I 1,...1,a k+1 ,...an (j) m.c. as a vector with an infinite number of components, indexed by non-positive integer tuples (a k+1 , . . . , a m ). Let d be the dimension of the vector space spanned by these s vectors (with meromorphic functions in D as coefficients). We then define the maximal cut as the linear basis of these s vectors, For the examples we considered in this paper, we find that d = n MI , the number of master integrals on the maximal cut. Furthermore, let S be the d×d matrix, whose element in the ith-row and jth-column is the i-th master integral evaluated on the subregion Ω b j of the maximal cut. We find explicitly that, for the examples we considered, S is the Wronskian matrix for the differential equation on the maximal cut.
It is also interesting to study the expansion of S near D = 4 (or any integer spacetime dimension). Define D = 4−2 . For example, the expansion is directly related to the -form of the differential equation [35,36] on the maximal cut level [43,49,50]. Suppose that for the i-th column of S, S i , this expansion reads where T i is the leading coefficient column vector, which is itself D independent. Let T = (T 1 , . . . , T d ) be the square matrix consisting of the leading coefficients. If the differential equation on the maximal cut reads, where A and B are independent, then by the -expansion of the differential equation. If T is invertible 1 then the new basis I = T −1 I satisfies the -form of the differential equation on the maximal cut, This is equivalent to the Magnus rotation in [38].

Maximal cuts without ISP
Our first example is the one-loop box integral in D dimensions with purely massless kinematics. Let k 1 , . . . , k 4 be the external momenta subject to the conditions k 2 i = 0 and 4 i=1 k i = 0. We define the two independent Mandelstam invariants by s = (k 1 + k 2 ) 2 1 If T is not invertible, then we may study the null vectors of T and the next-leading expansion coefficients of S to get the transformation matrix. and t = (k 1 + k 4 ) 2 . In order to simplify the problem we will study the maximal cut of the integral rather than the full integrated expression. There are no ISPs in this case, so we will instead consider integrals with arbitrary nonnegative powers a 1 , . . . , a 4 of the four propagators, The denominator factors are given by A constructive way of proceeding is to examine the Baikov representation. As discussed in the previous section, the D-dimensional scalar box can be written as the four-fold integral where F (z) is the Baikov polynomial. The maximal cut in arbitrary dimension is the quadruple cut realized by the replacement z −1 i → δ(z i ) for i = 1, 2, 3, 4. This cut localizes the box integral completely. The value of the maximal cut is thus basically determined by the Baikov kernel evaluated at the origin, The leading singularity evaluated in strictly integer dimensions has proven extremely useful when searching for and designing integrals that have uniform degree of transcendentality. Our compact analytic expression embodies the well-known result for the leading singularity in strictly four dimension, but also in odd dimensions, for example for D = 5, We can gain further insight by taking advantage of the maximal cut (3.4) to, for example, extract information about IBP relations for integrals with doubled propagators. It is straightforward to see that Upon comparison of eqs. (3.4) and (3.6) the common Gamma function can be dropped and therefore we immediately deduce the first very simple instance of an IBP relation, where integrals with fewer than four propagators are truncated. More generally, for generic values of the indices a 1 , . . . , a 4 ≥ 0, Therefore we can with almost no effort derive any desired IBP identity, for instance for integrals with tripled or several repeated propagators, simply by taking multiple derivatives and relating the resulting expression to eq. (3.4), I 2,2,1,1 = + (D − 6)(D − 5) st I 1,1,1,1 + · · · , (3.10) It is also worthwhile to investigate the dimensional dependence of the maximal cut, that is, examine the dimension shifting and dimensional reduction identities. The dimension shifting identity relates a D-dimensional integral with an extra-dimensional numerator insertion of µ 2 and a scalar integral in (D + 2) dimensions. At the level of maximal cuts we readily observe that The ratio of the maximal cuts in this equation is merely a polynomial function of D as it of course should be. In order to appreciate this fact, recall that the Gamma function satisfies the functional identity, Indeed, by iteration of eq. (3.13) we get the Pochhammer symbol, commonly referred to as the ascending factorial, 14) The dimensional reduction identity relating a (D + 2)-dimensional integral to a linear combination of D-dimensional scalar integrals is again the Gamma function property (3.13) in disguise. For example, Finally we look at differential equations and the maximal cut. Since we have access to the leading singularities we can make an educated choice for the normalization which leads to a canonical differential equation in the spirit of ref. [35]. Specializing to D = 4 − 2 and defining J = stI 1,1,1,1 we find the expected -form, where x = t/s. The exact same calculation goes through for D = 5 − 2 , where according to eq. (3.5) we should instead choose the master integral to be J = √ s √ t √ s + tI 1,1,1,1 . All results obtained in this section are consistent with the literature.

Maximal cuts with one ISP
The next simplest D-dimensional maximal cut is the integral with one ISP. In this section, we briefly review this case.
For example, we consider the D-dimensional two-loop three-point box-triangle diagram, with six propagators, see figure 1. The inverse propagators are and the external momenta satisfy k 2 1 = k 2 2 = 0, (k 1 + k 2 ) 2 = s. There are 2 × (3 − 1) + 2 × (2 + 1)/2 = 7 scalar products of the diagram, therefore the number of ISPs is 7 − 6 = 1. We may choose the ISP to be The Baikov variables are defined to be z i ≡ D i , i = 1, . . . , 7. The integrals under consideration are We may drop the argument D in the notation except for the discussion of dimension shift identities. We consider different mass configurations: • m 1 = m 2 = 0. In this case, the Baikov polynomial with maximal cut A generic integrand z k 7 /(z 1 . . . z 6 ), k ∈ N 0 , on the maximal cut reads, Note that so the integrand of (4.5) is a polynomial-valued total derivative. Hence The D-dimensional maximal cut vanishes and this implies that the massless boxtriangle integral is reducible (to integrals with fewer propagators).
The expression (4.16) can also be used for deriving the dimension shift identity, (4.23) which means the dimension shift identity on the maximal cut.
It is also interesting to study the differential equation on the maximal cut, where . . . stands for integrals with fewer propagators. Explicitly, we see that the Ddimensional cutĨ[1, 1, 1, 1, 1, 1, 0] m.c. solves the maximal cut part of this equation.
We can expand the D-dimensional maximal cut in (D = 4 − 2 ), where γ is the Euler-Mascheroni constant. From the leading coefficient in , we see that we may redefine the integral, such that the differential equation on the maximal cut has the form, where . . . stands for integrals with fewer propagators.
• m 1 = 0, m 2 = 0. Define x 1 = m 2 1 /s and x 2 = m 2 2 /s. Following similar steps as in the previous case, we definẽ . (4.29) We find that the D-dimensional maximal cut is . (4.30) Again, we can use it to study the IBPs, dimension shift identities, differential equations and -form, on the maximal cut.

Planar maximal cuts with two ISPs
Now we consider D-dimensional maximal cuts with two ISPs. In this case, the integral regions become less trivial and we can investigate each region seperately to get a complete set of D-dimensional maximal cut functions in kinematic variables.

Massless sunset
The massless sunset diagram. 2 All momenta flow to the right.
As a warm-up, consider the sunset diagram with massless internal propagators, see figure 2. We have the inverse propagators There are two ISPs, which we may choose to be where s = k 2 . Our object of interest is the integral Taking the Baikov variables to be z i = D i , equation (5.3) becomes where the integration region Ω is defined by F (z) ≥ 0. On the maximal cut the Baikov polynomial is which we can simplify by rescaling the Baikov variables to and our integral of interest on the maximal cut then reads , These regions are shown schematically in figure 3. The integration in each of the regions can now be carried out explicitly to give Note that the integrals only converge when certain restrictions on a, b, and D are satisfied, but we drop those as we are interested in the analytic continuation anyway. Assuming a and b to be non-negative integers, it follows from the reflection property of the gamma function, sin(πz)Γ(1 − z)Γ(z) = π, that So there is only one linearly independent function. (Here the coefficient field is set to be the meromorphic function field of D.) We thus set the D-dimensional maximal cut as the integration over the region I, (5.10) If we introduce the descending factorial 11) which is related to the Pochhammer symbol (3.14) by it is straightforward to see that .
In a similar way we readily find the dimension shift identity, (5.14) The differential equation for J[0, 0] m.c. is simple, 15) and is immediately in -form if we specialize to D = 2 − 2 .

Massless double box
The diagram is shown in figure 4. The inverse propagators are 16) and the external momenta satisfy k 2 1 = k 2 2 = k 2 4 = 0, (k 1 + k 2 ) 2 = s and (k 1 + k 4 ) 2 = t. We define χ = t/s. There are 2(4 − 1) + 2(2 + 1)/2 = 9 scalar products in loop momenta, hence we have two ISPs, The maximal cut of J[a, b] can be calculated by the integration of F (D−6)/2 z a 8 z b 9 . Consider the kinematic condition s > 0, χ > 0. The integration region defined by F ≥ 0 splits into four subregions: which are shown in figure 5 as the blue area. It is clear that on the four subregions the conditions µ 11 (z) ≥ 0 and µ 22 (z) ≥ 0 are satisfied.
The integrations over the first three subregions are straightforward, while the integration over the fourth subregion needs careful further splitting. The integration over Ω I where 2F1 is the regularized hypergeometric function, 2F1 (α, β, γ, z) = 2 F 1 (α, β, γ, z)/Γ(γ). This result is to be understood as an analytic continuation over D, a and b. Similarly the second integration over Ω II gives Although apparently this expression does not look symmetric in a and b, after a hyperge-ometric function transformation, the symmetry is manifest: The integration over Ω III gives m.c. .
These identities agree with the IBP output from FIRE [56][57][58][59], LiteRed [60,61], and Azurite [55]. We conclude that the function relations for J[a, b] m.c. provide the IBPs on the maximal cuts. Gauss' contiguous relations also imply dimension shift identities on the maximal cut. For example, for i =I, II, III, IV, Note that the cos(πD) factor in the region (5.24) does not affect dimension shift identities, since it is invariant under D → D + 2. We also remark that it is well-known that the solutions of recursive relations for the dimension shift identities may consist of hypergeometric functions [64]. We know that the two master integrals of the double box topology satisfy the differential equation on the maximal cut,

Double box with one massive leg
We can use the same method to consider the maximal cut of the double box with one massive external leg. The Feynman integral has the same inverse propagators (4.1), and the kinematic conditions are k 2 1 = m 2 1 , k 2 2 = k 2 3 = k 2 4 = 0. The two ISP are D 8 = (l 1 + k 4 ) 2 , D 9 = (l 2 + k 1 ) 2 and the diagram is shown in figure 6.
Define χ = t/s, κ = m 2 1 /s, x = z 8 /s and y = z 9 /s. The Baikov polynomial on the maximal cut reads Consider the kinematic region s > 0, χ > 0 and χ + 1 − κ > 0. Again, the integration region defined by F ≥ 0 splits into four subregions: while the integration over Ω II reads Again Gauss' contiguous relations provide IBPs on the maximal cut. We can also check that these maximal cut functions satisfy dimension shift identities.
Define the 2 × 2 matrix Explicitly, S is the fundamental solution matrix of the differential equations on the maximal cut, The leading coefficients of S in the limit D → 4 read Again, the transformation matrix is 44) and the new master integrals arẽ The differential equations turn into the -form,

Double box with two massive legs
The maximal cut of the double box with two non-adjacent external massive legs (k 2 1 = 0, k 2 3 = 0 or k 2 1 = 0, k 2 4 = 0) again consists of hypergeometric 2 F 1 functions, while the double box with adjacent two external massive legs (k 2 1 = 0, k 2 2 = 0) on the maximal cut yields Appell F1 functions. So in this subsection, we focus on the latter case.
The Feynman integral has the same inverse propagators (4.1), and the kinematic conditions k 2 1 = m 2 1 , k 2 2 = m 2 2 and k 2 3 = k 2 4 = 0. The two ISPs are D 8 = (l 1 +k 4 ) 2 , D 9 = (l 2 +k 1 ) 2 and the diagram is shown in figure 8. Define χ = t/s, κ 1 = m 2 1 /s, κ 2 = m 2 2 /s, x = z 8 /s and y = z 9 /s. The Baikov polynomial on the maximal cut is For example, we consider the kinematic regime where s > 0, κ 1 and κ 2 are small. Specifically, the latter condition means, It is then important to specify the sign of the expression inside square roots, so we pick up the condition (5.50).
The subregions are shown in figure 9. Explicitly, we can check that on the four subregions where F 1 is the Appell F1 function. The arguments w 1 and w 2 are defined as The integration over Ω II reads (5.56) The integration over Ω III is We can check that the last integration over Ω IV is dependent, m.c. . (i) m.c. , i =I, II, III, IV, satisfy IBP relations on the maximal cut. For example, where the 3 × 3 matrices are given in figure 10. Explicitly, for any i =I, II, III, IV, the (i) m.c. ) T solves the differential equations on the maximal cut. Furthermore, we find that the 3 × 3 matrix is the fundamental solution matrix for all these differential equations and the Wronskian is nonzero. The leading coefficients of S in the limit D → 4 are By a simple column operation, the transformation matrix can be chosen as So we can redefine the basis as The differential equations on the maximal cut forĨ is in the -form, where the matricesB χ ,B κ 1 andB κ 2 are given in figure 11. Again . . . stands for integrals with fewer propagators.
-κ +κ -κ -κ κ +κ Figure 11: Matrices for the -form differential equations on the maximal cut for the double box diagram with two adjacent massive legs.

Nonplanar maximal cuts with two ISPs
We are also interested in confirming that nonplanar integrals are equally amenable to the D-dimensional maximal cut procedure demonstrated for one-loop and planar two-loop integrals throughout the previous sections. More specifically, we analyze the maximal cut regions for the purely massless nonplanar double box, evaluate the two-fold integral in the post-hepta-cut degrees of freedom for each subregion separately, and finally investigate the properties of the resulting maximal cut expression.

Massless nonplanar double box
We adopt the conventions for the external kinematics from the preceeding section. Consider now a generic purely massless D-dimensional nonplanar double-box integral, specified by the seven inverse propagators, and two ISPs, which may be picked up as the conventional propagator-like forms, (Ω) where F (z 8 , z 9 ) = z 8 z 9 (s + z 8 ) (sχ − z 8 − z 9 ) 4s 2 χ(χ + 1) . (6.6) This expression for F parallels eq. (5.19). The domain of integration Ω is again defined as the region of the (z 8 , z 9 )-plane where the cut Baikov polynomial is nonnegative. As shown in figure 13, for the problem at hand, Ω is divided into five simple subregions. These subregions take the shape of triangles and rectangles, corresponding to the inequalities: We now calculate the double integral for each subregion separately. For example, the integration over region I is parametrized as follows, In this instance and also for the remaining four subregions the result is straightforwardly obtained in the form of 2 F 1 hypergeometric functions multiplied by kinematic invariants and some Gamma functions as expected. The result for region I can be written as

X[a, b]
Completely analogously, for region II, Next, we present the results for regions III and IV , Finally, the integration over region V yields a slightly more complicated expression involving multiple hypergeometric functions. As we shall see below, the fully simplified output from Mathematica can be further manually reduced to a very compact form. By inspection of eqs. (6.8)-(6.11) along with the expression obtained for region V , we see that the maximal cut of the nonplanar double box basically gives rise to Gauss' 2 F 1 hypergeometric functions with four distinct quadruples of indices, namely, Again, it is however trivial to reexpress two of the hypergeometric functions in terms of the remaining two. For instance, we can choose the integration over regions I and II as the principal results, and simply define the maximal cut function as An elementary manipulation of eq. (6.16) enables us to also include nonplanar double box integrals with doubled (or simply arbitrary powers of) propagators in the present setup. Here, however, we refrain for brevity from giving any examples along this direction. Instead, we verify that our maximal cut inherits the dimension shifting properties satisfied by the uncut integral. Explicitly, it can be shown that the maximal cut satisfies the raising recurrence relation, and also the lowering recurrence relation, These equations are true for any of the five subregions. Similar identities hold for all other nonplanar double box integrals. All integral relations inferred from the maximal cut are validated by FIRE [56][57][58][59] and Azurite [55]. Let us finally discuss differential equations obeyed by the nonplanar double box integrals and their relation to the maximal cut. From the IBP relations it has already been established that d dχ suppressing integrals with fewer than seven propagators. We have explicitly checked that our maximal cut indeed solves this differential equation, region by region. Moreover, the maximal cut functions (6.16) again form the Wronskian matrix S associated with this system of differential equations. Defining S as follows, the Wronskian determinant is found to be nonvanishing, m.c. are linearly independent functions of χ as previously anticipated, and furthermore confirms that the columns of S form the two fundamental solutions to eq. (6.25). The leading terms of S in the limit D → 4, again help us to identify a new set of master integrals in order to transform eq. (6.25) to -form. We may include a trivial redefinition of T and take the transformation matrix as which by eq. (2.12) implies that d dχ where the modified mastersJ 1 andJ 2 denote the following mixture of the FIRE basis integrals, The new differential equation is obviously in -form for D = 4 − 2 .
In summary, we have explicitly verified that all salient features of the maximal cut of the planar double box carry over immediately to the nonplanar double box. Our example here only covered the purely massless case, but we also have succesfully checked several configurations with massive external momenta.

Conclusion
In this paper we have presented a precise and consistent technique for evaluating generalized cuts of multiloop Feynman integrals, properly regularized in D dimensions. Our method relies on the Baikov representation, which makes the effect of taking these generalized cuts in arbitrary dimension manifest. We have given examples of the method for several integral topologies with various kinematic configurations, including the one-loop box, twoloop sunset, planar double box and nonplanar double box.
The simplest instance is the maximal cut of a k-propagator integral realized by the multivariate residue of the integrand at z 1 = · · · = z k = 0, the z i 's being the Baikov variables. In general the maximal cut leaves a multi-fold integral over a domain Ω defined by the intersection of the cut hyperplanes and the region A where the Baikov kernel F is nonnegative. The remaining integration may be carried out over any subregion of Ω with F = 0 on the boundary. We refer to the result as the maximal cut function on subregion j; this can be viewed as a generalization of the notion of a composite leading singularity.
The maximal cut function satisfies the same form of integral relations, such as IBPs and dimension shift identities, and differential equations as the uncut integral, region by region.
In our examples, the maximal cut functions are compact analytic expressions involving Gamma functions and (generalized-) hypergeometric functions. The integral identities on the maximal cut hence immediately correspond to relations among these special functions, namely recurrence relations and Gauss' contiguous relations.
One of the principal features of all our examples is that an integral topology with m master integrals has precisely m linearly independent maximal cut function series. For instance, the purely nonplanar double box, with real kinematics, gives rise to five subregions, but only two linearly independent maximal cut functions. This number agrees with the number of master integrals. We have explicitly shown that the linearly independent maximal cut functions form the Wronskian matrix S for the system of differential equations on the maximal cut. Moreover, we have in detail demonstrated that the leading terms of S provide a transformation matrix for differential equation into to canonical (epsilon) form near four dimensions.
From the viewpoint of the differential equation of Feynman integrals without cut, these maximal-cut functions form the fundamental solution set of the "homogenous" equation. To solve the complete differential equation system, can be understood as solving for an "inhomogenous" differential equation. Hence it is expected that these functions would appear in the complete expression of Feynman integrals.
This work brings inspiration for further advances along the direction of multiloop generalized unitarity with arbitrary spacetime dimension.
• In this paper, we simply consider the integration regions corresponding to the real loop momenta and find that it is enough to get the complete solutions for the differential equations on the maximal cut. For more complicated integrals, we may also consider complex loop momenta and integration regions for complex Baikov variables.
(For a region to be valid, we still require that on the boundary the Baikov polynomial vanishes, i.e., F = 0.) • It would be very interesting to study the maximal cut of elliptic Feynman integrals [42,[67][68][69][70][71][72][73][74][75] with arbitrary spacetime dimension. The 4D maximal cut of an elliptic double box was studied via Weierstrass elliptic functions [23]. The (D − 2)-iterative form of elliptic differential equations were studied in ref. [76][77][78]. We expect that our method combined with integrals over the fundamental domain of elliptic curves, would provide the maximal cut of elliptic Feynman integrals in a closed form of D.
• It is also interesting to see the non-maximal cut in D-dimension. We may either directly carry out Baikov integrals with non-maximal cut, or extend known maximalcut functions to the non-maximal cut results via the block-triangular differential equation.