On the maximal cut of Feynman integrals and the solution of their differential equations

The standard procedure for computing scalar multi-loop Feynman integrals consists in reducing them to a basis of so-called master integrals, derive differential equations in the external invariants satisfied by the latter and, finally, try to solve them as a Laurent series in $\epsilon = (4-d)/2$, where $d$ are the space-time dimensions. The differential equations are, in general, coupled and can be solved using Euler's variation of constants, provided that a set of homogeneous solutions is known. Given an arbitrary differential equation of order higher than one, there exist no general method for finding its homogeneous solutions. In this paper we show that the maximal cut of the integrals under consideration provides one set of homogeneous solutions, simplifying substantially the solution of the differential equations.


Introduction
The method of differential equations [1][2][3][4] is undoubtedly one of the most powerful techniques for the calculation of multi-loop Feynman integrals. The latter is based on the possibility of reducing a family of Feynman integrals to a small subset of basic integrals, the so-called master integrals (MIs), through the repeated use of integration by parts identities (IBPs) [5][6][7]. The IBPs themselves can then be used to show that the MIs fulfill systems of linear coupled differential equations in the external invariants, whose solution is usually much simpler than a direct integration over the loop momenta. For a review see [8,9].
We are normally interested in computing the master integrals as Laurent series in = (4 − d)/2. To this aim, there are two fundamental steps which must be achieved. First of all, given a sector with N coupled master integrals, one needs a way to determine the minimum number of coupled differential equations in the limit d → 4. This is important since it determines the effective degree of the differential equation that has to be iteratively solved. If all master integrals decouple (or if the homogeneous system takes a triangular form), we are left effectively with a series of first order differential equations in d = 4, which can be solved by quadrature. If, instead, some of the integrals remain coupled, one has to solve a higher order differential equation; naively one expects that the higher is the rank, the more involved the mathematical structure of the solutions will be. A way towards a systematical determination of the minimum degree of the coupled equations has been suggested in [10,11], where it was shown that the information useful to decouple some of the differential equations in d = 4 can be read off from the integration by parts identities close to d = 2 n, with n ∈ N.
The choice of master integrals is, of course, arbitrary and more recently it was shown that by properly selecting the basis of integrals one can simplify the form of the differential equations substantially, bringing them to a so-called canonical form [12]. One of the fundamental properties of a canonical form is that the homogeneous part of the system of differential equations becomes trivial in the limit → 0 where m is the vector of master integrals and x is a generic external invariant. This implies that the homogeneous solution for = 0 is a constant. Different algorithms have been proposed for the construction of a canonical basis, starting from the properties of the system of differential equations. Currently, they are all limited to special situations, such as a linear dependence on [13] or the dependence on a single kinematic variable [14][15][16].
Steps towards algorithms valid also in the case of several variables have been made in [17] and more recently and thoroughly in [18]. Both papers make use of Ansätze for the linear combinations of master integrals which fulfill canonical differential equations. While such approaches are often very useful, their applicability remains for now limited to those cases where all square roots can be removed and the alphabet can be completely rationalized. As it is well known, this is not always possible even when a canonical basis is known to exist, for a recent example see [20]. On the other hand, in [12] it was suggested that one of the crucial criteria for selecting candidate integrals which fulfill canonical differential equations is the study of their leading singularity, which can be done by inspecting their generalized cuts [9]. While a complete understanding of the concept of leading singularity and of a canonical basis are available only in the case of integrals that evaluate to multiple-polylogarithms and generalizations thereof 3 [19], it has been suggested that the study of the maximal cut should be the starting point to extend these considerations also to more complicated cases [20]. A crucial observation is that the maximal cut provides by construction exactly a solution of the homogeneous system of differential equations. The connection between unitarity cuts and differential equations is not new. In [21], for example, the second order differential equation satisfied by the two-loop massive sunrise graph was solved by inferring the homogeneous solution from the calculation of the imaginary part of the graph. More in general, this connection has been largely exploited in the so-called reverse unitarity [22] method, while a new way of solving IBPs using the information coming from the unitarity cuts has been proposed in [23]. Finally, similar conclusions to the ones drawn in this paper have been already exploited in the context of the DRA method based on dimensional recurrence relations [24].
A way to algorithmically find a homogeneous solution of a coupled system of equations is indeed very desirable. In fact, once we are confident to have reduced the degree of the differential equation to the minimum possible, the problem of finding an integral representation for the solution relies on the ability to solve its homogeneous part. Once a homogeneous solution is known, one can use Euler's method of variation of constants in order to build up the inhomogeneous solution. The goal of this paper is to show that the maximal cut can be used in many non-trivial examples as a powerful tool to compute the homogeneous solutions, in particular when the latter can be written in terms of complete elliptic integrals. One should recall, in fact, that given a coupled system of N differential equations, a complete solution requires finding N linearly independent solutions. In the case of elliptic integrals, one is usually left with systems of two coupled equations, which require therefore finding two independent solutions. While the maximal cut provides only one solution, once this is known, the second can be obtained by simply exploiting the properties of the complete elliptic integrals.
The rest of the paper is organized as follows. In Section 2 we set up the notation, summarize the general idea and give a prescription to compute the maximal cut of integrals with squared propagators. In Section 3 we apply the method to a simple one-loop example and show explicitly the one-to-one relation between the maximal cut and the solution of its homogeneous differential equation. This relation is true for any values of d and works of course both ways, which implies that it also provides us with a powerful tool to compute the maximal cut of any one-loop integral in d dimensions: this can be done simply by solving its homogeneous differential equations. In Section 4 we then move to more interesting two loop examples which evaluate to integrals over elliptic integrals. Finally we conclude in Section 5.

Maximal cut and differential equations
Let us consider a family of l-loop Feynman integrals with r different propagators and s irreducible scalar products, The integrals are functions of the dimensional regularization parameter d and of some external invariants which we call collectively x. First of all, we should define what we mean with cutting an integral and, in particular, with its maximal cut. It is well known that unitarity cuts can be used to study the discontinuity of Feynman graphs with respect to a given Mandelstam invariant. The original integral can then be reconstructed from its discontinuities by using dispersion relations [25][26][27]. For more recent developments in this directions see for example [28,29]. Here we will not enter into this fascinating and rich subject. The only thing we will need is an operative definition of maximal cut for an integral of the form (2.1). Cutting a propagator in a Feynman diagram means, loosely speaking, to force the particle propagating through it to be on-shell; mathematically, we can achieve this relatively easily if the propagator is raised to power one. In this case, we can simply substitute the propagator with a Dirac δ-function which forces the momentum of the propagator on-shell by imposing a constraint on the components of the loop momenta. We define the maximal cut of a diagram as the simultaneous cut of all its propagators and we will always assume such multiple cut to exist. In fact, if the set of on-shell conditions for a given diagram does not admit any solution, the diagram can be proven to be reducible, i.e. it can be trivially expressed as a linear combination of diagrams with a reduced number of propagators [30]. In the rest of our discussion will focus on diagrams corresponding to irreducible topologies, whose maximal cut is therefore a well defined object. As we will show in the following, we found it convenient to compute the maximal cut of two-loop graphs by first localizing one of the one-loop sub-diagrams, and then integrating directly the remaining δ-functions obtained cutting the second loop. By appropriately choosing which sub-loop to cut first, the computation can be substantially simplified. Moreover this procedure can easily be generalized for higher loops. We note here that an alternative way of performing the maximal cut is going through Baikov's representation for the loop integrals [31], as explained for example in [23].

Cutting squared propagators
It is useful to define what we mean with cutting a squared propagator. Indeed, at one loop one is always left with at most one master integral per topology, which can always be chosen with propagators raised at most to power one. On the other hand, in a generic multi-loop calculation, more integrals can remain independent and we might have to consider also integrals with squared internal propagators. In what follows, we will use two operative prescriptions to cut a squared propagator, which produce equivalent results. We will use these prescriptions to compute the maximal cut but, of course, they can be applied for the computation of any other cut.

The first prescription is based on integration by parts identities. Given any family of
Feynman integrals, we can use IBPs to write any integral with squared propagators in terms of similar integrals with propagators raised at most to unit powers and, possibly, scalar products at the numerator. Let us take such an integral I dot (d; x) and write where m j (d; x) are the n linearly independent master integrals which do not contain any squared propagator. The c j (d; x) are instead rational functions in the dimensions d and in the Mandelstam invariants x.
If we apply now a maximal cut on (2.2) we get where we used the fact that the maximal cut of the subtopologies is identically zero.
2. An alternative prescription to cut a squared propagator is as follows. Let us take an integral defined as where the propagators D k could in general depend on some masses m k , some or all of which could of course also be zero. We are integrating over l-loops and assuming that the integral has r propagators; the s-th propagator is squared. In order to compute the maximal cut of this integral, let us consider the associated integral If the limit is smooth, this provides us with a second operative prescription for computing the maximal cut of a graph with a squared propagator.

Differential equations and the maximal cut
Let us switch now to the connection between differential equations and the maximal cut defined above. Let us consider again a family of Feynman integrals like in Eq. (2.1), which is reduced to N independent master integrals m i ( ; x), ..., m N ( ; x). The master integrals satisfy a system of N coupled differential equations in the external invariants where the H ij (d; x) are the coefficients of the homogeneous system, while the G i (d; x) contain the dependence on the sub-topologies, which are simpler graphs with fewer propagators. Imagine now to perform a maximal cut on the master integrals m i (d; x), which corresponds to cutting all r propagators. As we discussed above, we only consider irreducible cases where all r propagators can be simultaneously cut. Clearly, if we imagine to apply our cutting procedure on Eq. (2.8) we will be left with where we use the notation Cut(m i (d; x)) for the maximal cut of the integral under consideration and we have obviously Eq. (2.10) should be obvious, since all integrals contained in G i (d; x) contain fewer propagators than the master integrals m i (d; x) and, consequently, applying the same cut on the latter must produce zero. We should stress here that, in general, the maximal cut exists only for complex values of the kinematical invariants. Eqs. (2.9, 2.10) are the central observations on which this paper is based. Similar conclusions in the context of the DRA method have been drawn in [24]. What they tell us is that the maximal cut of the master integrals m i (d; x) must satisfy the homogeneous part of the system (2.8). A remark here is in order. It is clear that if, for any reason, the maximal cut of any of the master integrals m i (d; x) turns out to be zero for a given value of the dimensions d, this implies that the master under consideration must decouple from the system (2.9) for this value of d. Indeed, while the trivial solution m i (d; x) = 0 is always a solution of the homogeneous system, the latter will in general admit other non-trivial solutions, which cannot be computed by evaluating the maximal cut. We will see an explicit example of this later on in Section 4.2.
Let us see how this works with a very simple example. We consider the one-loop massive bubble where p 2 = s. The latter satisfies very simple differential equations in s, m 2 1 and m 2 2 . The one in s, for example, reads where G(d; s, m 2 1 , m 2 2 ) is the inhomogeneous part which depends only on the tadpoles .
Eq.(2.13) defines also our integration measure, whose exact value will nevertheless be irrelevant in this context. The homogeneous part of the equation (2.12) reads (we use the suffix H to indicate that we are neglecting all inhomogeneous terms) and it admits the solution for generic d where c 1 is an irrelevant multiplicative constant. By the considerations above, Eq.(2.15) is also the d-dimensional maximal cut of the one-loop bubble. Clearly, in this case, the maximal cut coincides with the cut in the s-channel, which reads The cut is straightforward to compute in the frame p µ = ( √ s, 0) and one immediately obtains where c 2 is another constant whose exact value is irrelevant. Comparing Eq. (2.15) and Eq. (2.17) it is clear that, up to an irrelevant multiplicative constant, the solution of the homogeneous equation coincides with the maximal cut of the graph. A comment is in order. If a graph depends on more than one scale, like in the case above, the maximal cut provides the homogeneous solution of all differential equations in all Mandelstam invariants. In general, solving only the homogeneous equation in one of the invariants cannot capture the full dependence on all the remaining scales. For example, in the case of the one-loop bubble analyzed above, if we had solved only its homogeneous differential equations in ∂/∂m 2 1 , we would have not been able to determine the overall dependence on (−s) (1−d/2) in Eq. (2.17). Therefore, in order to get the full answer for the maximal cut we must solve all homogeneous differential equations in all Mandelstam invariants. While this example is straightforward, it contains most of the features of the more complicated examples that we will study below.
Having retained full dependence on d in the example above might look overly complicated. Indeed, for physical applications, we will be eventually interested in the solution of the system of differential equations as a Laurent series in = (4 − d)/2. The discussion above applies of course for any integer values of d as well and, in particular, for d = 4. If we go back to the general equation (2.8), put d = 4 − 2 and expand it in , we will quite in general end up with a system of equations in the form Since we are interested in solving (2.18) as Laurent series in , the crucial step consists now in solving the homogeneous system evaluated for = 0 Once the homogeneous solution is known, one can write order by order in an integral representation for the inhomogeneous solution, which will then be represented in terms of iterated integrals. Naively, we expect that solving (2.19) must be substantially simpler than retaining full dependence on d. While this is indeed true, there exists no general algorithm to solve (2.19) if more than one equation is coupled, even for d = 4. It is sometimes useful to rephrase (2.19) as a N -th order differential equation for any of the functions h j (x) and try to solve the latter. As of today, a limited number of two-loop examples are known which require the solution of second order differential equations. In all these cases, a solution could be found in terms of complete elliptic integrals of the first and second kind with complicated arguments and irrational prefactors. The best known example is the two-loop massive sunrise graph [10,21,[32][33][34][35], while more recently other, in some cases unrelated, examples have been worked out [20,[36][37][38][39]. In all these cases, once the homogeneous solution was determined it became possible to write useful integral representations for the complete result. Until now these results have been obtained with a case by case analysis. In this respect, very recently an interesting approach to solve a second order differential equation in terms of elliptic integrals has been proposed in [20]. The latter is based on the possibility of reparametrizing the differental equations in terms of a suitably chosen parameter, in terms of which its solution in the form of elliptic integrals becomes manifest.
We will follow here a complementary approach. We will show that, also in complicated cases which require the introduction of elliptic integrals, the solution of the homogeneous equation for = 0 can be rather simply obtained from the maximal cut of the integral evaluated in d = 4. We will work out explicitly different examples, including the double box considered in [20]. While this will allow us to produce equivalent results, we believe that our approach has an interest on its own, in particular since it can be at least in principle extended to any other arbitrarily complicated example. As we will see, the main limitation of our approach is the possibility of performing explicitly the integrals over the residual component of the loop momenta, after all δ-functions stemming from the propagators have been integrated out. This can potentially be a serious limitation. Nevertheless, we will show that in many cases this is not a problem and the calculation can be organized so to be left with only one residual one-dimensional integration. The complexity of this last integration depends, of course, on the integral under consideration. In simple cases which can be solved in terms of multiple polylogarithms, it can be usually performed in terms of simple rational functions or square roots of rational functions. Whenever the result involves elliptic function, instead, one can use a quite general change of variable in order to rephrase the result in standard form in terms of complete elliptic integrals. Finally, even if the last integrations cannot be performed analytically, this method allows to obtain relatively compact integral representations of one homogeneous solution of a complicated higher order differential equation, which is a very useful piece of information towards its complete solution.
In what follows we will see how the ideas described here work in practice in different examples. We will start analyzing more in detail a one-loop example, and move then to more interesting two-loop cases.

A one-loop three-point function
At one loop, for every topology, one finds always at most one master integral, which implies that one-loop integrals satisfy at most first order linear differential equations. It is then usually straightforward to integrate the homogeneous part of the differential equation by quadrature. According to the discussion above, this has to turn out to be identical to the computation of the maximal cut for the integral under consideration. Despite the simplicity of the calculations, we will do this here explicitly for an interesting example. This will allow us, on the one hand, to show some features of the procedure and, on the other, to obtain some useful formulas that we will recycle for the more interesting two-loop examples discussed below.
Let us consider a one-loop triangle with two massive internal propagators and three off-shell external legs where in total generality q 2 = q 2 1 = q 2 2 = 0. In this case it is very simple to compute the maximal cut of the graph, in particular in d = 4: We go to the reference frame q µ = ( q 2 , 0), q µ 1 = (E 1 , 0, 0, q 1z ), q µ 2 = (E 2 , 0, 0, −q 1z ) with and obtain immediately (we neglect at every step irrelevant overall numerical constants) where we definedk = | k| and ω a = k2 + m 2 a . Performing the integral in dk = ω a /k dω a we get where we have fixedω a = (q 2 + m 2 a − m 2 b )/2/ q 2 . Given that there exists a real or complex value for the momenta where the δ-function in (3.4) has support 4 , we are left with Now, following our argument, this cut must be proportional to the solution of the homogeneous differential equations satisfied by Tri(4; q 2 , q 2 1 , q 2 2 , m 2 a , m 2 b ). It is very simple to derive the differential equations satisfied by the latter in all external invariants. As exemplification we write down here the homogeneous differential equation in q 2 for generic d, which reads (again, we use the suffix H to indicate that we are neglecting all inhomogeneous terms in the differential equation) where C(q 2 , q 2 1 , q 2 2 , m 2 a , m 2 b ) is a cumbersome rational function. In the simpler case where m b = m a the latter reads . (3.7) Putting d = 4 in Eq. (3.6), one sees that the entire dependence on the masses m a and m b disappears and the solution reads where c is an integration constant. This is of course in perfect agreement with the result of the maximal cut calculation, Eq. (3.6). Two comments are in order here. First of all, we could have of course solved Eq. (3.6) retaining full dependence in d. The result is relatively compact and reads Again, up to an overall constant, this is also the result of the calculation of the d-dimensional maximal cut of this triangle. A second interesting feature of our result is that, in the limit d → 4, it does not depend either on m a or m b . This has some important consequences. Let us consider a similar integral, but with one of the two massive propagators squared, for example such that Tri 2a (d; q 2 , q 2 1 , q 2 2 , m 2 a , m 2 b ) = ∂ ∂ m 2 a Tri(d; q 2 , q 2 1 , q 2 2 , m 2 a , m 2 b ) .
Following our prescriptions for computing the maximal cut of this integral, we define simply  The very same thing is true if we put a dot on the other massive propagator. We should stress here once more that this result is valid only for d = 4 identically. It is clear from Eq. (3.9) that taking a derivative w.r.t. m 2 a or m 2 b would produce an overall factor (d − 4) which goes to zero as d → 4.
A similar conclusion can be drawn inspecting the integration by parts identities. For ease of writing, we consider again the case m b = m a . Using IBPs one can show that which, given Eq.

Two-loop examples
At two-loop order the situation is, indeed, much more interesting. In this case, one often has to consider topologies of Feynman integrals which are reduced to two or more master integrals. The latter will therefore satisfy systems of coupled differential equations, whose solution becomes much less straightforward. In this case, the first step consists in finding a solution of the homogeneous system. Once the latter is known, one can proceed using Euler's variation of the constants to write down an integral representation for the full solution. In what follows we will considered two examples of two-loop Feynman graphs which fulfill second order differential equations and show explicitly how the maximal cut allows us to find the required homogeneous solutions.

A non-planar two-loop triangle
Let us start considering a two-loop non-planar triangle with two off-shell legs and four massive propagators which enters the non-planar corrections to H+j production through a massive fermion loop. We define our integral family as follows q p 1 p 2 = I a 1 ,a 2 ,a 3 ,a 4 ,a 5 ,a 6 ,a 7 a 7 <0 with p 2 1 = 0, p 2 2 = 0 and q 2 = (p 1 + p 2 ) 2 = s = 0 . The integrals (4.1) can be reduced to two master integrals that we choose as Neglecting all subtopologies, the latter satisfy the following two sets of coupled differential equations in the variables x = −s/m 2 and y = −p 2 2 /m 2 , d dx d dy where B x,y (x, y) and D x,y (x, y) are 2 × 2 matrices that do not depend on , given by (4.4) In d = 4, ( = 0) the systems become d dx and they can be rephrased as second order differential equations for one of the two master integrals. For I 1 they read and Instead of trying to solve these equations directly, we compute the maximal cut of I 1 in order to determine a first homogeneous solution. I 1 can be written as (4.8) The integral in D d l is a one-loop box with four massive propagators and three off-shell external legs , (4.9) where p 1 and p 2 are defined above, while q 3 = k − p 1 and q 4 = −k − p 2 . In order to compute the maximal cut of I 1 , it is convenient to start by first cutting the one-loop box Eq. (4.9). Since we have already spelled out the relation between maximal cut and homogeneous differential equation, we can avoid to compute this cut by direct integration over the δ-functions. Instead, we derive the homogeneous differential equations satisfied by the latter in all external invariants and solve them for d = 4. The result provides us with the maximal cut, up to an overall irrelevant constant. We obtain where we have defined t = (p 1 − q 3 ) 2 and u = (p 2 − q 3 ) 2 . Substituting this into (4.8) and localizing the contour for the remaining two propagators we are left with , (4.11) where we used the condition q 2 3 = q 2 4 = 0. In order to compute the second integral, we go to the center of mass frame and perform the integration in k 0 andk and obtain where we used the definition of the complete elliptic integral of the first kind , for x ∈ C and (x) < 1 , .
In terms of the adimensional variables x and y, the maximal cut becomes (neglecting once more overall constant factors) (4.14) By direct computation one can check that F 1 (x, y) solves both second order differential equations (4.6) and (4.7). Finally, the properties of elliptic functions allow us to write a second independent solutions as (4. 15) We observe that in the limit p 2 2 → 0 we have which provides the homogeneous solution for the differential equations of the corresponding non-planar two-loop triangle with only one off-shell leg. This solution can be then analytically continued to all physically relevant phase-space region and can be used to write down a one-fold integral representation for the two-loop triangle. This problem will be considered elsewhere [40].

An elliptic planar double box
As a last non-trivial example, we consider the planar double-box computed in [20], which corresponds to the integral family = I a 1 a 2 a 3 a 4 a 5 a 6 a 7 a 8 a 9 a 7, 17) and the propagators are defined as We stress that, since we are interested in the six-denominator topology only, we restrict the indices a 7, 8,9 to negative values. The external kinematics is chosen to be p 2 1 = p 2 2 = p 2 3 = 0, p 2 4 = (p 1 + p 2 + p 3 ) 2 = m 2 h and the Mandelstam invariants are given by The IBPs reduction of the integral family (4.17) returns four master integrals, which we choose to be The master integrals I i fulfill systems of first order differential equations in the kinematic invariants. In addition, it can be verified that in d = 4 the differential equations for I 4 are completely decoupled. Therefore, in the four-dimensional limit we can restrict ourselves to study the homogeneous systems for the first three master integrals, which read where x ∈ x = {s, t, m 2 h , m 2 } and a ij (x) are rational functions in the Mandelstam invariants. Following our argument, we expect that computing the maximal cut of the three master integrals we will get, by construction, one set of homogeneous solutions of the system (4.21).
Interestingly enough, it is easy to show that cutting all propagators produces a further simplification of the system. In fact, we observe that I 3 contains as sub-loop the one-loop triangle with one dotted massive propagator studied in Section 3. There we saw that the latter has vanishing maximal cut in d = 4. Hence, as a direct consequence of Eq.(3.12), we immediately get The vanishing of Cut (I 3 ) implies that, when the systems (4.21) are evaluated on the maximal cut, the last equation drops out and we are left with two coupled homogenous first order differential equations for I 1 and I 2 , Cut (I 2 (x)) = a 21 (x) Cut (I 1 (x)) + a 22 (x) Cut (I 2 (x)) .

(4.23)
A comment on the relation between the solutions of the homogenous systems (4.21) and the existence of a vanishing maximal cut is in order. Since the differential equation for I 3 is decoupled from I 1 and I 2 , a zero maximal cut provides an obvious solution to it. However, Cut (I 3 ) = 0 does not prevent the differential equation to admit other non-trivial solutions, which are not captured by the maximal cut. Nevertheless, this is not a problem since the differential equation is decoupled and can be solved independently by quadrature. We would like to stress here that the decoupling of the differential equation for a master integral in four dimensions and the vanishing of its maximal cut are indeed closely related. In fact, precisely a vanishing maximal cut could be seen as the hint of the decoupling of the differential equation for the corresponding master integral, since only a decoupled equation would be automatically satisfied by a zero solution without imposing strong constraints on the maximal cuts of the other two master integrals.
The systems of equations (4.23) can then be rephrased as a second order differential equations for the maximal cut of one of two master integrals. For instance, if we choose I 1 and we differentiate with respect to the internal mass, we have and other three similar equations can be determined by taking derivatives with respect to s, t and m 2 h .
One of the two independent solutions of this set of second order differential equations can be found by direct computation of the maximal cut of I 1 . As in the previous example, we start by computing the maximal cut of one of the two sub-loops from the solution of the corresponding differential equations in d = 4 and then integrate over the second loop momentum by solving explicitly the constraints imposed by the remaining δ-functions. In this respect we observe that, in principle, one is completely free to choose the order in which the loop momenta are localized. Nevertheless, a wise choice can substantially simplify the remaining integrals. In this particular case, the two sub-loops correspond to a box with two adjacent off-shell legs and a triangle with three off-shell legs. If we started by cutting the box, we would be left with an integral over the four variables parametrizing the second loop momentum and only two δ-functions constraining them. Therefore, in order to obtain a one-fold integral representation of the solution, we would need to perform an additional non-trivial integration. On the other hand, if we localize the triangle integral first, the action of the three remaining δ-functions would directly provide an expression of the maximal cut in terms of a one-dimensional integral.
We start then from the cut of the internal one-loop triangle which, using the results of Section 3, reads .

(4.25)
Therefore, the maximal cut of I 1 can be written as where we have used that, when the remaining propagators are cut, (k − p 1 − p 2 ) 2 = m 2 . Before discussing the evaluation of the last integral, we would like to make a technical remark. When the number of on-shell conditions to be imposed is sufficiently large, the kinematic restrictions of Minkowski space could suggest that no non-trivial solution of the maximal cut exists. Therefore, one needs to adopt a practical prescription in order to obtain a result in some kinematic region and then continue it to the region of interest through a consistent relaxation of the assumptions used in the intermediate steps of the calculation. For this specific case, we found that an effective prescription consists in assuming a negative internal mass, m 2 ≤ 0 . Under this assumption, we evaluate the integral over the loop momentum k in the rest frame of the massive external particle p µ 4 = (m h , 0) , by parametrizing the massless momenta p µ 1 and p µ 2 as The energies of the two particles and the relative angle between their three-dimensional momenta are expressed in terms of the Mandelstam invariants as . (4.28) If we go to polar coordinates and we decompose the loop momentum as k µ = (k 0 ,k sin θ cos ϕ,k sin θ sin ϕ,k cos θ), (4.29) where we have again defined | k| =k, we can easily check that the three additional cutconditions imply . (4.31) The δ-function δ(k 2 −m 2 ) can be then used to perform the integral overk by fixing its value tō k = k 2 0 − m 2 . In this way, the constraint imposed by δ(2k · p 1 ) reduces to a linear equation in cos θ, which yields to . (4.32) 18 Finally, we evaluate the integral over the angle ϕ by using the constraint imposed by the last δ-function, which is satisfied by We observe that, in order for these solutions to range in 0 ≤ ϕ ± ≤ 2π, we must restrict the integration over k 0 within the region with In this way, we obtain a one-fold integral representation of the maximal cut of I 1 of the type , (4.36) where we have introduced two additional roots defined as The integral (4.36) can be cast into the canonical form of a complete elliptic integral of the first kind through the standard change of variables which yields to a solution of the second order differential equation (4.24) of the form with (4.40) We observe that F 1 has a smooth behavior in the limit of vanishing internal mass m 2 → 0, which reproduces the correct result for the maximal cut of a six-denominator box with massless propagators. Finally, using the properties of elliptic integrals, it is simple to find a second independent solution of Eq. (4.39) by changing the argument ω of the elliptic function to 1−ω, We verified explicitly that Eq. (4.39, 4.42) satisfy the second order differential equation (4.24), and of course also the corresponding ones in the other Mandelstam invariants s, t and m 2 h . We can obtain an alternative (but equivalent) representation of the solutions as follows. We recall that the complete elliptic integral of the first kind is defined as , for x ∈ C and (x) < 1 , (4.43) and it fulfills a second order differential equation in the form Since any second order differential equations admits only two independent homogeneous solutions, one must have that where c i , i = 1, 2, 3, 4 are numerical (possibly complex) constants. Of course, since K(x) develops a branch cut when x > 1, one should assign a small imaginary part to x, which determines the sign of the imaginary parts of the coefficients c i . We find, for example, for 0 < x < 1 and x → x + i δ, that the following relations are satisfied This means that we can equally well use either the F 1 and F 2 defined above or the two new solutions Finally, it is interesting to compare our result to those of [20]. There, the solution was found suitably redefining the Mandelstam variables in terms of an additional dimensionless parameter α, with respect to which the solution of the equations became straightforward. For α = 1 one should recover the standard solution. Indeed for α = 1, the argument of the elliptic integrals found in [20] reduces precisely to 1/ω, as in Eqs. (4.47).

20
Differential equations are an invaluable tool for computing multi-loop Feynman integrals. For any practical purposes, the complexity of a system of differential equations is largely dictated by the number of coupled differential equations that must be solved in the limit d → 4. First of all, one needs to find a complete set of homogeneous solutions; once they are known, Euler's method of the variation of constants can be used to reconstruct the complete inhomogeneous solution. As a matter of fact, as soon as we are left with a system of two or more coupled irreducible differential equations, there exists in general no algorithm to find a complete set of homogeneous solutions; the possibility of solving the system depends therefore on a caseby-case analysis.
Building upon ideas developed in the context of the study of the differential equations satisfied by the two-loop massive sunrise graph [21], we showed that the maximal cut of any given (irreducible) Feynman integral provides us precisely with one set of homogeneous solutions of the differential equations satisfied by the latter. This one-to-one relation can of course be inverted and, in some case, one might want to use the solution of the homogeneous differential equations as tool to compute the maximal cut of a graph. This is particularly useful at one-loop, since every one-loop integral satisfies at most a first order differential equation, whose homogeneous part can always be solved by quadrature. In this sense, there is no need to compute the maximal cut of any one-loop integral by direct integration over the loop momenta, as the latter can always be obtained by solving its homogeneous differential equation, for any arbitrary values of the dimensions d.
For higher loops this is not true and we are left in many cases with higher order equations. In this case, it turns out to be much easier to compute the maximal cut then solving the homogeneous equation directly. The maximal cut provides us with one set of homogeneous solutions. When dealing with second order differential equations, the second solution can then be obtained in closed analytic form with standard methods. Far from being a simple curiosity, we showed explicitly that this observation can be turned into a very powerful tool for solving second (and possibly higher) order differential equations. In different non-trivial two-loop examples, in fact, we showed how the maximal cut can be very easily computed in d = 4 space-time dimensions, providing therefore a formidable piece of information towards the solution of coupled systems of differential equations. A crucial step in this direction is recognizing that the calculation of the maximal cut of a multi-loop integral can be suitably split in the calculation of the maximal cuts of the relevant sub-loops. A wise choice in the order in which sub-loops are cut can simplify substantially the calculation of the maximal cut.
Finally, we believe the results of this paper are interesting also in view of the extension of the concept of canonical basis to Feynman integrals that do not evaluate to iterated integrals over d-log forms only. Suppose, in fact, to have a N × N system of differential equations in the form which is one of the fundamental properties a canonical basis should fulfill. A thorough investigation of these aspects is beyond the aim of this paper and will be subject of study in the near future.