Integration by parts identities in integer numbers of dimensions. A criterion for decoupling systems of differential equations

Integration by parts identities (IBPs) can be used to express large numbers of apparently different d-dimensional Feynman Integrals in terms of a small subset of so-called master integrals (MIs). Using the IBPs one can moreover show that the MIs fulfil linear systems of coupled differential equations in the external invariants. With the increase in number of loops and external legs, one is left in general with an increasing number of MIs and consequently also with an increasing number of coupled differential equations, which can turn out to be very difficult to solve. In this paper we show how studying the IBPs in fixed integer numbers of dimension d=n with $n \in \mathbb{N}$ one can extract the information useful to determine a new basis of MIs, whose differential equations decouple as $d \to n$ and can therefore be more easily solved as Laurent expansion in (d-n).


Introduction
Dimensionally regularised [1][2][3] Feynman integrals fulfil different identities, among which the most notable ones are the so called integration by parts identities (IBPs) [4,5]. Given a family of Feynman integrals, the IBPs can be used to write down a large system of linear equations with rational coefficients and which contain the Feynman integrals of that family as unknown. By solving algebraically the system a large number of apparently different Feynman integrals can be expressed in terms of a much smaller basis of independent integrals dubbed master integrals (MIs). In realistic application the number of such equations can grow very fast, requiring the use of computer algebra in order to handle the complexity of the resulting expressions. There are different public and private implementations which allow to perform the reduction to MIs in a completely automated way [6][7][8][9] based on the so-called Laporta algorithm [10,11].
The IBPs can be used to prove that dimensionally regularised Feynman integrals fulfil linear systems of first order differential equations in the external invariants [12][13][14][15][16]. Considering a Feynman graph which is reduced to N independent MIs, by direct use of the IBPs one can derive a system of N coupled linear first order differential equations for the latter, which can be rephrased as and N -th order differential equation for any of the MIs. Supplemented with N independent boundary conditions, the system of differential equations contains all information needed for a numerical or analytical calculations of the MIs. Indeed, in the general case, the analytical solution of an N -th order differential equation is a very non trivial mathematical problem.
It has been observed that, in many cases of practical interest, a substantial simplification of the problem occurs when studying the behaviour of the system of differential equations as the space-time dimension parameter d approaches 4, which is also the physically relevant case. Usually we are indeed not interested in an exact solution for the MIs as functions of d, but instead in the coefficients of their Laurent expansion for d ≈ 4. In [17,18], and in many subsequent applications of the differential equation method, it was shown that it is often possible to choose a basis of MIs such that the differential equations take a simpler triangular form in the limit d → 4. If this is possible, the problem of integrating the system of differential equations simplifies substantially, reducing de facto, at every order in (d − 4), to N subsequent integrations by quadrature. Experience showed that, whenever such a form is achievable, the differential equations can be integrated in terms of a particular class of special functions, the multiple polylogarithms (MPLs) [17,19,20]. The latter have been studied extensively by both mathematicians and physicists and routines for their fast and precise numerical evaluation are available since some time [21][22][23]. Disclosing their algebraic properties allowed furthermore the development of very powerful tools for the analytical manipulations of these functions [24][25][26].
More recently it has been shown [27,28] that in many of these cases a basis of MIs can be found, such that the system of differential equations takes a particularly simple form, commonly referred to as canonical form. The system is said to be in canonical form if the regularisation parameter d can be completely factorised from the kinematics, appearing as an explicit (d − 4) factor in front of the matrix of the system. In addition, the coefficients of the matrix must be total-differentials of logarithms of functions of the external invariants (i.e. they are said to be in d-log form). A canonical basis is particularly convenient as it allows a straightforward integration as series expansion in (d − 4) and, due to the d-log form of the coefficients, it integrates directly to MPLs of uniform transcendental weight. Criteria for the construction of candidate canonical integrals have been presented in [28] and developed in detail, for example, in [29]. In the special cases in which the differential equations depend only linearly on the dimensions d, the Magnus algorithm can be used to perform a rotation of the system to a canonical form [30]. For a recent application of the algorithm see [31]. A completely different approach based on Moser algorithm [32] has been developed for the univariate case in [33] and discussed also independently in [34]. Another interesting approach is based on the properties of higher order differential equations fulfilled by the individual master integrals [35]. In spite of all this impressive progress, a fully automated algorithm, working also in the multivariate case, is still missing. It has nevertheless become clear that, if one can find a basis of MIs whose differential equations become triangular as d → 4, it is often (but not always) relatively easy to bring it in canonical form by removing the undesired terms in the differential equations [36].
Indeed, different examples are known where neither finding a canonical basis nor a triangular one is possible. In all these cases the master integrals cannot be integrated in terms of simple multiple poly-logarithms only [37][38][39][40]. It becomes therefore a very interesting problem that of finding a set of criteria to determine whether, given a Feynman graph, there exists a basis of MIs whose system of differential equations becomes triangular in the limit for d → 4. This would provide a way to easily classify all possible diagrams where this is not possible and that would therefore be expected to evaluate to more complicated classes of function. Aim of this paper is to show that a large amount of information about the possibility of achieving such decoupling can be extracted by studying the IBPs for fixed even numbers of dimensions, i.e. in the limit d → 2 n, where n is any integer number. Tarasov-Lee shifting identities [41,42] allow, in fact, to directly relate the structure of the differential equations in any even number of dimensions with the physical case d = 4. Indeed, since Feynman integrals are usually divergent at d = 2 n, this limiting procedure must be carried out as a Laurent series in (d − 2 n).
The rest of the paper is organised as follows. In Section 2 we review the use of integration by parts identities for the reduction to master integrals and we summarise the main results of the differential equations method. In Section 3 we outline the central idea of the paper. We show in particular what kind of information can be extracted by solving the IBPs as Laurent series for d → 2 n and how to use it to simplify the system of differential equations. In Section 4 we then apply these ideas explicitly to many different examples of increasing complexity. Some comments on the method are given in Section 4.7, where we try as well to point out some relevant open issues. Finally we conclude in Section 5. In Appendix A we compare our method with the Schouten identities introduced in [39], while in Appendix B we show how to shift a system of differential equations of an even number of dimensions.

Integration by Parts identities and Differential Equations
Let us consider an l-loop scalar Feynman integral depending on P independent external momenta p i where k i are the loop momenta, D i = q 2 i − m 2 i are the propagators and S i = k j · p l are irreducible scalar products. In view of the discussion below, we will usually denote the integrals of a topology as in (2.1), keeping explicitly only the dependence on the dimensions d and on the powers of denominators and scalar products. In dimensional regularisation, for every integral of the form (2.1) there always exists a value of the space-time dimensions d, such that the integral is convergent 2 . Necessary condition for the convergence of an integral is the integrand be zero at the boundaries. This condition can be mathematically rephrased as where the v µ m are any of the external or internal momenta v µ m = {k 1 , ..., k l ; p 1 , ..., p P } . Such an identity is called an integration by part identity or IBP. It is clear that in this way l(l + P − 1) IBPs can be established for each integrand. Upon explicitly evaluating the derivatives and contracting with the momenta v µ m , new integrals belonging to the same topology (i.e. integrals with the same set of denominators) are generated. In particular, each IBP identity can relate integrals with (s − 1), s and (s + 1) powers of scalar products, and (t + r) or (t + r + 1) powers of propagators. Notice that, by contracting with v µ m , new reducible scalar products can be generated, which could then simplify some of the denominators producing integrals belonging to any of the (t − 1) sub-topologies of the original graph.
It has been shown [10,11,16] that the system of IBPs, which appears in general to be over-constrained, can instead be solved allowing to express most of the integrals as linear combinations of a small subset of basic integrals, dubbed master integrals (MIs). Indeed, as for any algebraic basis, the choice is not unique and, by suitable changing the basis, one can substantially simplify the calculation of the integrals, as we will discuss later in this paper.
Let us consider now a Feynman graph (or a topology) characterised by a set of propagators D i and irreducible scalar products S i . All integrals will depend on the space-time dimensions d and on the external invariants x ij = p i · p j , where p i are as usual the external momenta. Let us assume, for definiteness, that by solving the system of IBPs 3 all integrals for the given graph can be expressed in terms of a basis of N independent MIs I i (d; x ij ) with i = 1, ..., N , which are of the same form of (2.1). Differentiating with respect to any of the external invariants x ij amounts to differentiating with respect to linear combinations of the external momenta p µ i [14]. Therefore, by acting with these differential operators directly on the integrands of (2.1), one produces linear combinations of integrals belonging to the same Feynman graph and to its sub-topologies. The latter can again be reduced to MIs, generating in this way a system of N linear first order differential equations with rational coefficients in any of the invariants x ij . Suppressing the dependence on the sub-topologies, which can be considered as a known inhomogeneous term in a bottom-up approach, the homogeneous part of the system can always be written as where the coefficients c ij (d; p 2 ) are simple rational functions of the dimensions d and of the external invariants x ij . We can rewrite the system in matrix form as where we introduced the vector of master integrals I(d; x ij ) and the matrix of the coefficients A(d; x ij ).
The system (2.3) is, in the general case, coupled and can therefore be rephrased as an N -th order differential equation for any of the MIs I i (d; x ij ). In most practical applications we are interested in determining the MIs as Laurent expansion for d ≈ 4 (2.5) By expanding both left-and right-hand-side of (2.3) one is left with a chained system of N differential equations where, at any order α, the previous orders can only appear as inhomogeneous terms.

An optimal basis of master integrals
It has been shown by Tarasov and Lee [41,42] that the value of a Feynman integral in d space-time dimensions can be directly related to that of the same integral in d − 2 or d + 2 space-time dimensions. This implies that, if all MIs of a given graph are known as Laurent expansion in any even number of dimensions, d = 2 n, then the coefficients of their series expansions in d = 4, I i (4; x ij ) in (2.5), can be obtained as linear combinations of the I (α) i (2 n; x ij ). For more details see for example [37,39] and the discussion in Appendix B. Indeed, changing the basis of MIs changes the form of the matrix A(d; x ij ) in equation (2.4). An interesting problem is therefore how to define an optimal basis of MIs in order to simplify as much as possible the system (2.3). Since we are interested in computing the MIs as Laurent expansion in (d − 4) (or, in general, in (d − 2 n)), an obvious simplification would occur if we could decouple some of the differential equations, at least in the considered limit. In particular, given a system of N coupled equations, one could think of classifying the complexity of the latter by determining the minimum number of differential equations that cannot be decoupled in the limit d → 4 (or more generally d → 2 n).
At this point it is useful clarify more precisely what we mean with decoupling in this context. Let us consider a 2 × 2 coupled system of differential equations 4 where I(d; x) = (I 1 (d; x), I 2 (d; x)) is the 2-vector of unknown functions, A(d, x) is a 2 × 2 matrix, d are the space-time dimensions and x is a variable the functions depend on 5 . Assume now for simplicity that the functions I 1 (d; x), I 2 (d; x) are finite in the limit d → 4 and that the matrix A(d; x) does not contain any explicit poles in 1/(d − 4). Assume finally that, in the limit d → 4, the matrix A(d; x) has non-zero non-diagonal entries, and therefore that the system is coupled in the limit d → 4. Of course, by expanding the entries of A(d; x) in (d − 4) we can write our system as It is now clear that if, by any means, we can find two independent solutions to the 2 × 2 system , v 2 (x)) and (w 1 (x), w 2 (x)), then we can define the new vector J (d; x) through the rotation such that the differential equations satisfied by J (d; x) assume the form i.e. they become trivial in the limit d → 4. The matrix G(x) can be of course arbitrarily complicated, as it contains the solutions of a second order differential equation. In this case we would have of course achieved a decoupling, but at the price of having to solve a coupled system of differential equations, which is in the general case not possible. On the other hand, what we are really interested in is to determine whether a basis of MIs exists such that some of the non-diagonal terms of the matrix A(d; x) become zero in the limit d → 4, and such that this basis can still be reached from our starting basis only through IBPs (i.e. without having to solve a coupled system of differential equations!). What this means in practice is that, if such a basis existed, then the rotation matrix G would assume a very simple form, namely it would contain only rational functions of the external invariants x ij (and of the dimensions d). This new basis would therefore fulfil a system of differential equations where some (or all) of the MIs decouple in the limit d → 4, and still it would be a system of linear differential equations with rational coefficients only. In this respect we note that, for all known cases of MIs which can be integrated in terms of multiple polylogarithms, a change of basis in the sense described above can be found and the system of differential equations can be put in triangular form as d → 4 where T (4; x ij ) is a triangular matrix and does not depend on the dimensions d. From the point of view of the classification outlined above this corresponds to the easiest case, where all equations decouple in the 4 Again, we neglect the inhomogeneous terms everywhere. 5 In the case of Feynman integrals x represents a generic Mandelstam variable.
limit d → 4 and, effectively, the problem reduces to many independent integrations by quadrature. Finding a basis in this form is often a first step towards a canonical basis in the sense introduced in [28]. We want to stress here that of course all these considerations apply in the very same way for any integer number of dimensions d → n (even or odd). Unfortunately a change of basis of this kind cannot always be found. Several cases are known where the system cannot be completely triangularised and instead at least two differential equations remain coupled. In these cases MPLs turn out not to be enough for describing the solution and the class of functions must be enlarged to include also elliptic generalisations of the latter [43]. It is unclear whether this will be the end of the story, since cases where 3 or more coupled equations survive are relatively easy to find, as we will show in the following. What still appears to be missing is a criterion to determine, given a Feynman graph, what is the minimum number of equations which cannot be decoupled. Together with simplifying as much as possible the problem at hand, this could also give a hint to which class of functions are required for describing the solution.

Reading the IBPs in fixed numbers of dimensions
In order to find a possible working criterion to determine the minimum number of coupled differential equations we should go back to think how the differential equations are derived. We saw that differentiating a master integral with respect to the external invariants produces new integrals belonging to the same Feynman graph. By using the IBPs one can then reduce these integrals to MIs, ending up with a system of differential equations. If we start with N master integrals we will obtain in general a coupled system of N linear differential equations. The fact that the N differential equations are coupled can be seen, in this respect, as due to the linear independence of the N master integrals in d dimensions.
As we already discussed, for any physical application we are interested in computing Feynman integrals as Laurent series in (d − 4) or, more generally, in (d − 2 n) with n ∈ N. Of course, different integrals have different degrees of divergence, i.e. their Laurent expansion starts at different orders in (d − 2 n). For any value of the dimensions, nevertheless, the maximal divergence can be computed in dimensional regularisation and depends only on the topology of the graph under consideration (i.e. on the number of loops, of external legs etc.). We can therefore imagine to first generate the IBPs in d dimensions and then expand them as a Laurent series in (d − 2 n), obtaining in this way a chained set of systems of IBPs, one for every order in (d − 2 n). It is clear that, by construction, at every order in (d − 2 n), the homogeneous part of each system will be identical, while the inhomogeneous part will contain the previous orders of the expansion (and the sub-topologies, that we will neglect throughout). If we limit ourselves to the first order of the expansion, i.e. the one corresponding to the highest pole in (d − 2 n), the system of equations that we are left with is equivalent to the original system of IBPs where d is fixed to be d = 2 n, and corresponds to the sole homogenous system. Now, it is very well known that upon fixing the number of space-time dimensions in the IBPs to an integer value it may happen that some of the equations degenerate and, in particular, that some of the integrals that used to be linearly independent for generic values of d, become linearly dependent from each other. From the point of view of the differential equations satisfied by the master integrals, if some of the integrals were to become linearly dependent in the limit d → 2 n, one would expect that those masters should not bring any new information in that limit and it should therefore be possible to decouple them from the system of differential equations as d → 2 n. Let us try to state this point more precisely. As exemplification we consider a topology that is reduced to 2 master integrals which we call I 1 (d; x) and I 2 (d; x), where d are the dimensions and x is a generic Mandelstam variable. Neglecting the sub-topologies the system of differential equations that they satisfy can be written as where the c ij (d; x) are rational functions. Let us now follow the argument above and generate the IBPs fixing d = 2 n. Let us assume that, by solving this simplified system, one of the two master integrals becomes linearly dependent from the other one and the new IBPs produce the relation I 2 (2 n; x) = b(x) I 1 (2 n; x) , (3.2) where b(x) is a rational function of the Mandelstam variables 6 . Equation (3.2) implies that in d = 2 n one of the two master integrals becomes linearly dependent in the sense of the IBPs. According to the argument above we would therefore expect to be able to decouple the two differential equations in this limit. In order to see this it is useful to ask ourselves how such a relation can emerge from the original d-dimensional IBPs. Let us imagine that upon solving the IBPs for generic d, we can find a d-dimensional relation expressing a given integral of the graph under consideration, say K(d; x), as a linear combination of the two masters and such that It is clear that, if this is the case, the IBPs which would generate this identity for generic d, would instead generate (3.2) once d is fixed to be d = 2 n. These relations are precisely what we are looking for. To refer to the latter we will often use throughout the paper the notation where it should be understood that, in general, this does not mean that the combination above is really of order O(d − 2 n), but simply that it becomes zero upon setting d = 2 n in the IBPs. Note that, of course, can only produce corrections of order O(d − 2 n) due to (3.3). We will see many examples of these relations in the sections below. Naively, the fact that only one integral is linearly independent for d = 2 n would require that the integral itself should satisfy a first order differential equation as d → 2 n. Finding a basis of master integrals for which this first order equation emerges is equivalent to finding a basis which decouples the system (3.1). To this aim let us perform the following rotation of the master integral basis The system (3.1) under this rotation becomes Equations (3.7) do not look particularly illuminating at first glance. We claim nevertheless that these equations are precisely what we were looking for. The basis J 1 (d; x), J 2 (d; x) defined in (3.6), in fact, has been chosen in order to exploit the linear dependence of the two master integrals in the limit d → 2 n. In this limit the IBPs tell us that J 1 (d; x) is by construction suppressed by a factor (d−2 n) and therefore decouples from the problem. We expect therefore that the differential equation for the latter should decouple in this limit or, in other words, that If this is true then upon expanding the system of differential equations as Laurent series in (d − 2 n) one can, at every oder, first solve the differential equation for J 1 (d; x) by quadrature, and then use this as an input for the second equation. A rigorous mathematical proof of equation (3.8) is outside the scope of this paper and we will limit ourselves to show explicitly how this works in practice with several examples of different complexity. The considerations above can be easily generalised to N master integrals I 1 ,..., I N . In this case one starts with a system of N coupled differential equations. By solving the IBPs for d = 2 n one can then verify how many of the master integrals become linearly dependent in this limit. Assuming that N − M integrals remain independent, this means that M relations like (3.3) can be found, say 9) and the b ij (d; x) are as always rational functions of the dimensions and of the Mandelstam variables 7 . As for the previous example we will often write these relations as where once more we imply that these combinations become zero upon setting d = 2 n in the IBPs. As before we define b ij (x) = lim d→2 n b ij (d; x) and, following the same reasoning, we can then try to rotate the basis of master integrals to Under the rotation (3.11), we expect the M integrals J 1 (d; x), ..., J M (d; x) to decouple from the remaining independent integrals in the limit d → 2 n, as in (3.7). One must be cautious here on what is intended by decoupling. According to the arguments above, upon the change of basis (3.11), we expect the system of differential equations to split into two blocks in the limit d → 2 n, one M × M and the other (N − M ) × (N − M ). This would correspond, order by order in (d − 2 n), to an M -th plus an (N − M )-th order differential equation, unless for some other reason internally the two blocks of differential equations further decouple in this limit. On the other hand, for the Feynman graphs that we considered so far (see for example Sections 4.3 and 4.6), even a stronger claim can be made. In these cases the rotation (3.11) not only splits the system into two blocks, as described above, but it also produces an explicit (d − 2 n) in front of the whole M × M block originated from relations (3.10). This explicit overall factor allows to effectively reduce the problem to the solution of one single (N − M )-th differential equation, plus M integrations by quadrature. The reason for this behaviour is still partly unclear and deserves further study. Summarising, the discussion above brings us to the following conclusion. Given a topology with N master integrals which fulfil a set of N coupled differential equations in d space-time dimensions, the study of the IBPs in fixed numbers of dimensions, say d = 2 n, provides a tool to determine how many master integrals can be decoupled from the differential equations as d → 2 n. Of course, the arguments given above are partly oversimplified and we have not provided here any rigorous mathematical proof. The structure of the differential equations can be in general very involved and, instead of embarking on complicated mathematical arguments, we prefer to show explicitly how this ideas can be simply applied to different cases of increasing complexity. In the next section we will start off considering simple examples where, by fixing the number of dimensions to an even integer value, only one master integral remains linearly independent and therefore the problem can be reduced to the solution of one linear differential equation. We will then move to more interesting cases where, even in fixed numbers of dimensions, more than one master integral remain linearly independent and one cannot avoid the problem of solving higher order differential equations which give rise to more complicated mathematical structures.

Explicit examples
In the previous section we outlined the main ideas behind this paper. We have discussed that the IBPs might degenerate in the limit of fixed (even) integer numbers of dimensions d → 2 n, such that some of the master integrals become effectively linearly dependent from each other. While this is a very well know fact, we argued that this degeneracy, if present, can be used in order to simplify the system of differential equations satisfied by the master integrals. In this section we will present many explicit examples of this very simple idea. We will start studying the two-loop sunrise graph with one massive and two massless propagators, Section 4.1, and a two-loop triangle with three off-shell legs, Section 4.2.
In both examples there are only two master integrals and by studying the IBPs in fixed even numbers of dimensions, one relation can be found, allowing to decouple the differential equations in that limit. We will then consider the case of the two-loop massive sunrise, Section 4.3, and of a non-planar two-loop triangle, Section 4.4. In both cases not all equations can be decoupled, and a minimal bulk of two differential equations remains coupled, giving rise to elliptic functions. We will then study the case of a two-loop massive triangle with three master integrals, Section 4.5. Here, similarly to the non-planar two-loop triangle, there are three master integrals. In this case, nevertheless, the differential equations can be completely decoupled and the solution can be written in terms of MPLs. Finally, as a last example, we will move to the three-loop massive banana graph 4.6. In this case, we will study different possible mass-arrangements of increasing complexity, showing how the number of master integrals changes consequently, and how our method allows to determine easily which subset of master integrals can be immediately decoupled from the differential equations.

The two-loop sunrise with one massive propagator
Let us start off by considering the case of the two-loop sunrise with one massive and two massless propagators. We define the following set of integrals belonging to its Feynman graph where p 2 = s is the momentum transfer. The integration measure is defined as and the explicit form of the function C(d) is not relevant for the considerations below. Note that this Feynman graph does not contain any sub-topology. We keep explicit only the dependence on the spacetime dimensions d and on the powers of the denominators and scalar products, which will be important for what follows. Performing a usual reduction through IBPs one finds two independent MIs, which can be chosen as I 1 (d; s) = I(d; 1, 1, 1, 0, 0) , Using the methods outlined in the previous sections we can now derive the differential equations fulfilled by I 1 and I 2 in the momentum squared s. This step can be performed automatically using, for example, Reduze 2 [9], and we end up with the following 2 × 2 linear system As one can immediately see, the equations are coupled for any even value of the dimensions d 8 .
It is well known that these integrals can be computed as a series expansion in d → 4 in terms of HPLs only, see for example [44]. Let us then try and use the ideas outlined in Section 3 in order to decouple the differential equations in the limit d → 4. First of all note that, in the limit d → 4, both master integrals are U V divergent and in particular they both develop a double pole Moreover it is easy to show that any integral of the form (4.1) can at most develop a double pole in (d − 4). Equipped of these consideration, let us now produce the IBPs for this Feynman graph as described above but, instead of solving them keeping the full dependence from the parameter d, we can set d = 4. 9 As we discussed in detail in the previous section, this is equivalent to expanding the IBPs in Laurent series, and considering the first of the chained systems of equations obtained, namely the one corresponding to the double pole in (d − 4). Following the arguments of the previous section, we would expect to find a degeneracy of the two master integrals in d = 4, which should then allow us to decouple the two differential equations. As expected the two masters (4.3) become linearly dependent In the notation of Section 3 we can write this relation as  .7) can be seen as a real relation between the highest poles of the two master integrals. This relation, which is naturally derived from the IBPs only, can be easily verified by computing the highest poles of the two master integrals. A very simple exercise gives , (4.9) in agreement with eq. (4.7). The overall normalisation of (4.9) is of course arbitrary and it has to do with the choice for the integration measure (4.2). Let us now try and exploit this relation in order to simplify the system of differential equations (4.4). We perform the change of basis from the "standard" MIs I 1 (d; s), I 2 (d; s), to the new MIs defined as Note that in this case the first of the two masters in (4.10) has only a single pole in (d−4) due to exactness of relation (4.7). As second master integral we can choose any of the two and here we performed simply a random choice picking I 1 (d; s). Choosing I 2 (d; s) would indeed bring to equivalent results. Deriving the differential equations for the new basis we find Equations (4.11) confirm the discussion in Section 3 and can therefore be seen as of the main result of this paper. Let us have a closer look at these two equations and compare them to (4.4). We note immediately that the equations are not in canonical form. On the other hand, the matrix of the system does become triangular in the limit of d → 4, and in particular the master J 2 appears in the the differential equation for the integral J 1 multiplied by an explicit factor (d − 4), as predicted in (3.8). For any practical purposes this is enough, since it means that one can expand the differential equations as Laurent series in (d − 4) and, order by order, first solve the differential equation for J 1 by simple quadrature, and then use this result as input for the differential equation for J 2 , which can in turn be solved by quadrature. Needless to say, this procedure can in principle be iterated up to any order in (d − 4). A comment is in order. In this simple example the relation found by studying the IBPs in d = 4 can be interpreted as an actual relation between the double poles of the two master integrals. Very often the first poles of randomly chosen MIs are either constants or very simple rational functions. One might therefore naively think that, by evaluating explicitly the poles of a given set of MIs, one could simply look for simple relations among the latter. Such relations would indeed contain only simple rational functions. It is well known, though, that in several cases the poles of a master integrals can be represented entirely through its sub-topologies. If this is the case, a relation between the poles of the masters would be useless, as it would not bring any new information as far as the master integrals are concerned. In order to achieve a decoupling one must therefore use a relation which is contained in the IBPs and as such represents an effective degeneracy of the master integrals in the limit d → 2 n, with n ∈ N.

Simplification of the differential equations in d = 2
It is interesting to see what happens by repeating the same analysis for the graph (4.1) in d = 2 instead of d = 4. Again, an easy analysis of the two master integrals (4.3) shows that they both develop a double pole in d = 2, which in this case is of IR origin As before, one can easily see that also in this case all integrals of the form (4.1) can develop at most a double pole in (d − 2). We can proceed and generate the IBPs for generic d and then expand them as Laurent series, this time in (d − 2), starting from 1/(d − 2) 2 . We are then left with a series of chained systems of IBPs, each for a different order in (d − 2). As in the previous case, we can now focus on solving the first system, corresponding to the double pole. This again is equivalent to considering the original system of IBPs and simply fixing d = 2. Upon doing this one immediately sees that once more the two MIs become linearly dependent (4.14) Proceeding as above, we can choose as new basis of MIs  Deriving the differential equations for J 1 and J 2 one finds immediately (4.16) Again, as expected from our analysis in d = 4, we see that differential equation for J 1 decouples from the one for J 2 in the limit d → 2, respecting the same pattern described in equation (3.8). Once more for every practical purposes this is enough to reduce the solution of the system of differential equations to iterated integrations by quadrature.

A two-loop triangle with three legs off-shell
Let us consider now a massless two-loop three-point function with three legs off-shell. The problem has been widely studied in the literature, mainly in the context of vector boson pair production [45][46][47][48], and it is well known that this Feynman graph can be reduced to two independent MIs, which can be integrated in terms of MPLs only. We define the Feynman graph as follows I(d; n 1 ,n 2 , n 3 , n 4 , n 5 , n 6 , where p 2 1 = m 2 1 , p 2 2 = m 2 2 and q 2 = (p 1 + p 2 ) 2 = s. We used Reduze 2 in order to reduce this graph to two independent MIs We can then proceed and derive the differential equations for these two MIs. As always we neglect the sub-topologies throughout. In this particular case the latter are simple two-loop corrections to massless two point functions which have been known analytically for a very long time.
The homogeneous part of the differential equations in the momentum transfer s reads where we defined the polynomial The equations are coupled in the limit d → 4. Again, as for the sunrise studied in Section 4.1, the integrals can develop at most a double pole in (d − 4). Instead of performing a complete Laurent expansion of the IBPs, we generate them and then fix explicitly d = 4 before solving them. This is enough to check whether the two MIs degenerate in this limit. By solving the IBPs one finds that this is precisely the case and the following relation is extracted where again we used the notation introduced above, indicating that the combination ( Note, nevertheless, that this time the relation is not exact, differently from (4.7), since the sub-topologies might in general contribute modifying (4.22). Eq (4.22) is anyway sufficient to decouple the homogeneous part of the differential equations. We proceed as above and define the new basis  Deriving the differential equations satisfied by (4.24) we find Again, as expected, we see that the differential equation for J 1 contains an explicit factor (d − 4) multiplying the second integral J 2 . The result is consistent with the general structure described in Section 3. Once more we want to stress that, as expected, in this case all MIs can be integrated in terms of MPLs only.

The two-loop massive sunrise
In the previous sections we considered two simple examples of 2 × 2 systems where the two equations could be decoupled in the limit d → 4, such that the problem could always be reduced to integrations by quadrature. As we already discussed this is not always possible and the first known case where at least two differential equations remain coupled is the two-loop massive sunrise. The two-loop massive sunrise graph is defined as follows I(d; n 1 , n 2 , n 3 , n 4 , This integral has been studied widely in the literature and in particular a lot of attention has been devoted to the differential equations that it fulfils. In the general case where all three masses assume different values, a normal reduction through IBPs shows that all integrals can be expressed as linear combinations of 4 independent MIs, which can be chosen to be In [15] it was shown that these integrals fulfil a coupled system of 4 linear first order differential equations in d dimensions. The system remains coupled in the limits d → 2 n, where n ∈ N. It was lately shown in [49], using algebraic geometry methods (and as such a priori orthogonal to the IBPs) that the scalar integral I 1 (d; s) satisfies a second-order Picard-Fuchs differential equation in d = 2. This suggested the possibility of finding a proper change of basis of MIs, in the sense of the IBPs, such that two of the four differential equations satisfied by the latter would decouple in the limit d → 2. Since the four MIs in (4.27) are finite in d = 2, it appeared natural to try and obtain the decoupling of the differential equations by finding new relations among the MIs, valid strictly only for d = 2. Such relations can be found using the so-called Schouten Identities introduced in [39]. In that reference the Schouten identities are introduced and the case of the two-loop sunrise with different masses is worked out in detail. It is shown that, as expected, in d = 2 only two master integrals are linearly independent. This allowed to recover the second order differential equation found in reference [49] in a completely independent manner. In this section we will show that those results can be even more easily re-obtained using the methods described in this paper, and namely by solving the IBPs for the massive sunrise in d = 2. The Schouten identities can be imagined as a tool for extracting this piece of information from the IBPs and are, in this respect, equivalent to the study of the IBPs in fixed number of dimensions. We will show another example of this equivalence in Appendix A.
Since the algebra in this case is rather heavy due to the large number of scales, we will only report the result of the solution of the IBPs in d = 2, referring to [39] for their use to simplify the system of differential equations. As already discussed above, solving the system with d = 2 is in general easier and, as expected, two of the four MIs degenerate, leaving only two linearly independent MIs. By choosing as MIs I 1 (2; s) and I 2 (2; s), we find the following additional relations (as everywhere else we neglect the sub-topologies for simplicity) As we discussed in Section 3, we expect such relations to come by d-dimensional IBPs with an overall factor 1/(d − 2). Indeed by studying the reduction to MIs in d dimensions it is easy to find the following two relations Relations (4.29) can be compared with the corresponding formulas (3.14) and (3.15) of [39]. It is easy to see that they are identical in the limit d → 2, the only difference being the absence of the terms coming from the sub-topologies, which we are neglecting here. These two relations (and in particular their limiting value as d → 2) can be used, as described in Section 3, in order to decouple two of the four differential equations of the two-loop massive sunrise graph, by choosing as new basis of master integrals We do not give the explicit form of the differential equations, referring to [39] for further details. In comparing, note that the basis presented here differs from the one in [39] by the absence of sub-topologies and by orders O(d − 2). Furthermore we want to stress, in relation to the discussion in Section 3, that using this basis produces an overall factor (d−2) in front of the two differential equations for J 3 (d; x) and J 4 (d; x). This implies that one has, at every order in (d − 2), only one second order differential equation (needed to solve the block of J 1 (d; x) and J 2 (d; x)), plus two integrations by quadrature (required to determine J 3 (d; x) and J 4 (d; x)).

A two-loop non-planar crossed vertex
As a further application, let us consider a two-loop non-planar crossed vertex with two massive propagators. This graph is topologically completely unrelated to the two-loop sunrise and was studied thoroughly in [50].
There it was shown that it can be reduced to three MIs, which would therefore be expected to satisfy a system of three coupled differential equations. In [50] a basis of MIs was found such that one of the three differential equations decouples from the other two in the limit d → 4. This reduced effectively the problem to that of solving, for every order in (d − 4), a second order differential equation, plus an integration by quadrature for the third MI. In this section we would like to study this Feynman graph with our method and show that the decoupling found in [50] comes as well from a degeneracy of the master integrals in d = 4 which can be read off directly from the IBPs. Following [50] we define the Feynman graph as follows I(d; n 1 , n 2 , n 3 , n 4 , n 5 , n 6 , n 7 ) = m m with p 2 1 = p 2 2 = 0 and (p 1 + p 2 ) 2 = s. It is easy to verify that this topology can be reduced to 3 MIs, for example Let us derive the differential equations in the momentum transfer s, neglecting as everywhere else all sub-topologies. We get Looking at these equations we see immediately that I 3 is already decoupled from the other two in the limit d → 4. This means that, with this randomly chosen basis, the problem is reduced to that of solving, at every order in (d − 4), a coupled system for I 1 and I 2 . With the explicit solution for the latter, one can then obtain I 3 integrating its differential equation by quadrature. It would be then interesting to know whether this decoupling is also due to a degeneracy of the MIs in d = 4. Moreover it would be even more interesting to verify whether a new basis could be found, for which the differential equations become completely triangular as d → 4, reducing even further the complexity of the problem.
Following the recipe described above, we can try and solve the IBPs for this Feynman graph for d = 4. A word of caution is required here. The three MIs selected above have Laurent expansions in (d − 4) which start at different orders, in particular one can easily find (for example using sector decomposition [51]) that the first two masters are finite, while the third develops a cubic pole (4.34) Nevertheless this poses no practical obstacle to the applicability of the method presented in this paper. As we already discussed in general, by expanding the system of IBPs in Laurent series around d = 4, we will get, at every order in (d − 4), an independent system of differential equations whose homogeneous part (i.e. the one containing the order of the MIs under consideration) has always the same form, while the non-homogenous part will of course change and, in particular, depend on the previous orders of the expansion (and on the sub-topologies, which we neglect). What we are interested in is, indeed, only the homogeneous part of this system. Fixing d = 4 is therefore enough in order to determine whether, for any order of the expansion, the MIs become linearly dependent. Upon doing this we find only one relation between the three masters which reads where, as always, we mean that this combination becomes zero when we set d = 4 in the IBPs. Equivalently, one can also proceed in a more formal way, expanding all IBPs in Laurent series starting from the triple pole up to the finite piece, and supplementing the piece of information on the highest poles of the MIs (4.34). Upon doing this, one obtains four chained systems of IBPs (one for every oder in (d− 4)), which can be solved bottom-up starting from the one corresponding to the highest pole. Since the first two masters are finite, the first three systems give no information on the latter, while the fourth system (corresponding to the finite piece of the MIs) produces the relation     1, 1, 1, 1, 1, 1, 1). If we had solved the IBPs in d dimensions, this integral would have been of course expressed in terms of the three masters (4.33). Solving the system in d = 4 instead does not allow to express this integral in terms of the other three, but this comes with no surprise and can be very well understood in terms of the degeneracy of the system of IBPs in this limit. By studying explicitly the integral I(d; 1, 1, 1, 1, 1, 1, 1) it is easy to see that it is also finite in d = 4, namely I(d → 4; 1, 1, 1, 1, 1, 1, 1) = O(1) , −→ I (−1) (4; 1, 1, 1, 1, 1, 1, 1) = 0 .
With this piece of information one recovers again relation (4.35), which was found by simply solving the system of IBPs in d = 4. Since only one relation has been found, which moreover involves all three masters I 1 , I 2 , and I 3 , we have no way to decouple more than one MIs from the system. What we mean here is that, since in system (4.33) only the differential equations for I 1 and I 2 are coupled, if we had found one relation but involving only I 1 and I 2 , we could have used it to decouple this block of the system. An example of this is given in Section 4.5. This is not the case and we therefore expect the minimal number of coupled integrals in d = 4 to be two, giving rise to a second order differential equation, as for the case of the two-loop massive sunrise, see Section 4.3.
As an exercise, we can try to change basis also in this case exploiting the piece of information found in (4.35). We expect to end up with a new system of differential equations, where nevertheless again two out of three equations are coupled as d → 4 (and as such practically equivalent to (4.33)), showing that the system can not be further simplified. Let us introduce the new basis Deriving the differential equations and neglecting all sub-topologies we get which is again a system of three differential equation, two of which remain coupled in the limit d → 4, giving rise to a second order differential equation for one of the two coupled masters. 10 We note that the new system (4.38), compared with the previous one (4.33), has a slightly different structure. As for the previous cases that we analysed, once we switch to the new basis defined through the IBPs degeneracy (4.35), the differential equation for the new master J 3 develops an explicit factor (d − 4) in front of the other two masters J 1 and J 2 , as predicted in equation (3.8).

A two-loop massive triangle with three master integrals
Before moving to a three-loop example, let us try and see what happens in a case similar to the one studied above, i.e. a Feynman graph reduced to three master integrals, but where the system of differential equations can be completely triangularised as d → 4. Let us consider the following two-loop massive triangle I(d; n 1 ,n 2 , n 3 , n 4 , n 5 , n 6 , n 7 ) = ✲ with two legs off-shell, namely P 2 = (p 1 + p 2 ) = s, p 2 1 = 0, p 2 2 = q 2 . This graph has been studied in the context of the QCD corrections to H → Zγ in [52,53]. Similarly to our previous example it is reduced to three master integrals, such that we are dealing with a system of three differential equations. We start from a randomly chosen basis of master integrals 1, 1,1, 1, 0, 0, 0) , I 2 (d; s, q 2 ) = I(d; 1, 1, 2, 1, 0, 0, 0) 1, 1, 1, 2, 0, 0, 0) . (4.40) The masters depend on three variables s, q 2 and m 2 , and therefore on two independent ratios. For simplicity we will consider only the differential equations in s, while all considerations done here work identically for the differential equations in the other variables. In order to simplify as much as possible the formulas we write explicitly only the order zero of the homogeneous differential equations in (d − 4), which is also the bulk which we need to simplify. The equations read where the dependence from q 2 is left as implicit in the master integrals for ease of notation. One can immediately see that only two of the differential equations are coupled. One should in principle first solve the 2 × 2 coupled system for I 2 (d; s) and I 3 (d; s), and then, with the latter as an input, one could attempt to solve the differential equation for I 1 (d; s) by quadrature. Let us try now and study the IBPs in the limit d → 4. By solving them as discussed in the previous sections one sees that the master integrals I 2 (d; s), I 3 (d; s) become linearly dependent in this limit and one finds the relation As for the case of the non-planar triangle studied in Section 4.4, we find only one relation, while we have three master integrals. One of the three masters nevertheless is already decoupled, and moreover relation (4.42) involves only I 2 (d; s, q 2 ) and I 3 (d; s, q 2 ), which are precisely the two coupled integrals. We expect therefore this to be enough to decouple the system. We define the new basis where the 1/m 4 has beed added for dimensional reasons. Deriving the differential equations for this new basis, and keeping again only the first order in (d − 4) we find As expected the system of differential equations becomes triangular and, in particular, the equation for the new integral J 3 (d; s), defined through relation (4.42), decouples from J 2 (d; s), following the usual pattern of equation (3.8). One can then proceed, order by order in (d − 4), integrating by quadrature first the differential equation for J 3 (d; s), then the one for J 2 (d; s) and finally the one for J 1 (d; s). As a last comment we want to stress that, if we had considered the system in ∂/∂q 2 , the same change of basis (4.43) would have indeed been sufficient to triangularise that one as well.

The three-loop massive banana graph
As last example let us consider a more complicated three-loop graph. We choose the three-loop massive banana graph, which is the natural three-loop generalisation of the two-loop massive sunrise. In the most general case this Feynman graph depends on the momentum squared p 2 = s and on four different masses m 1 , m 2 , m 3 and m 4 I 4 (d; n 1 , n 2 , n 3 , n 4 ,n 5 , n 6 , n 7 , n 8 , n 9 ) = ✲ ✫✪ ✬✩ p where the subscript 4 indicates that the four masses are all different. In the two-loop case there are 4 MIs when the 3 masses have all different values, which in turn degenerate to 2 MIs in the case of equal masses. On the other hand we saw that, irrespective of the explicit values of the internal masses, one is always left with only two independent MIs in d = 2 (4.28). This allowed us to decouple two of the four MIs from the differential equations in the limit d → 2 and prove that the scalar amplitude satisfies a second order differential equation in the this limit.
It would therefore be interesting to verify whether a similar behaviour can also be seen in the three-loop banana graph. Since in the general case with four different masses the algebra becomes very cumbersome, we will consider different cases of increasing complexity, namely increasing at every step the number of different internal masses and check how many MIs are found in d dimensions and how many can be decoupled in the limit d → 2.

The equal-mass case
Let us start considering the equal-mass case. We use the following notation I 1 (d; n 1 , n 2 , n 3 , n 4 , n 5 , n 6 , n 7 , n 8 , n 9 ) = I 4 (d; n 1 , n 2 , n 3 , n 4 , n 5 , n 6 , n 7 , n 8 , n 9 ) m4=m3=m2=m1=m , where the subscript "1" indicates now that all masses have the same value. Running a reduction to MIs with a code of choice it is easy to check that there are 3 independent MIs in d dimensions which can be chosen to be I 1 (d; s) = I 1 (d; 1, 1, 1,1, 0, 0, 0, 0, 0) , I 2 (d; s) = I 1 (d; 2, 1, 1, 1, 0, 0, 0, 0, 0) , d; 3, 1, 1, 1, 0, 0, 0, 0, 0) . (4.46) The differential equations in the momentum transfer for these three MIs read and we can easily check that, in spite of the fact that I 3 does not appear in the first equation, the system is still coupled as d → 2. Trying to solve the system of IBPs in d = 2 shows no further degeneracy and therefore we conclude that the system cannot be further simplified with our method. Having a system of three coupled first order equation means that we can rephrase it as a third-order differential equation for any of the three masters, and in particular for the scalar amplitude I 1 (d; s). The fact that the scalar amplitude fulfils a third-order differential equation is in agreement with the findings in [49]. Deriving the third-order differential equation satisfied by I 1 (d; s) we find where the d-dimensional third order differential operator D and all sub-topologies are neglected as always. In the limit d → 2 the differential operator simplifies to which is in agreement with [49].

The case of two different masses
Let us move now to a slightly more general case and let the masses take two different values. There are two possible arrangements, which we call I A 2 and I B 2 , defined as follows I A 2 (d; n 1 , n 2 , n 3 , n 4 , n 5 , n 6 , n 7 , n 8 , n 9 ) = I 4 (d; n 1 , n 2 , n 3 , n 4 , n 5 , n 6 , n 7 , n 8 , n 9 ) m3=m2=m1=ma,m4=m b , I B 2 (d; n 1 , n 2 , n 3 , n 4 , n 5 , n 6 , n 7 , n 8 , n 9 ) = I 4 (d; n 1 , n 2 , n 3 , n 4 , n 5 , n 6 , n 7 , n 8 , n 9 ) The two configurations are intrinsically different and it makes sense to look at the two cases separately.
A) In configuration A a reduction to MIs for generic d gives 5 independent MIs which can be chosen as 1, 1, 1,1, 0, 0, 0, 0, 0) , I A 2 (d; s) = I A 2 (d; 2, 1, 1, 1, 0, 0, 0, 0, 0) , 3, 1, 1, 1, 0, 0, 0, 0, 0) , 2, 2, 1, 1, 0, 0, 0, 0, 0) . A natural question at this point would be how many MIs degenerate in the two cases in the limit d → 2, and therefore what is the order of the differential equation satisfied by the scalar amplitudes I A 1 (d; s) and I B 2 (d; s) respectively. A naive expectation, based on the two-loop sunrise, would be to see in cases A and B, 2 and 3 MIs decouple respectively, such that the problem would reduce to the solution of a third-order differential equation, as in the equal-mass case. Unfortunately this naive expectation is not satisfied and we find that, by solving the IBPs in d = 2, in both cases four MIs remain independent, corresponding in principle to a fourth-order differential equation for the scalar amplitude in both mass-configurations. On the other hand, it is interesting to see that in both configurations, in spite of the different number of MIs in d dimensions, the problem can be reduced to an equation of the same order (i.e. four) in d = 2.
Neglecting the sub-topologies we find in configuration A the following relation which allows to express the fifth master integral in terms of the previous four In configuration B, instead, there are two different relations, which can be used to two express I B 5 and I B 6 in terms of the other four MIs in d = 2. We do not report the explicit solution of the IBPs in d = 2 which looks rather cumbersome. As for the case of the two-loop sunrise, these identities are originated from d-dimensional IBPs which degenerate in the limit d → 2 due to an overall factor 1/(d − 2). There are many of these relations, but only two of them are linearly independent in the limit d → 2, and they read (keeping only the first order in (d − 2)) We stress again that relations (4.53), (4.54) and (4.55) are not exact since all sub-topologies have been neglected throughout. These relations can be nevertheless used in order to derive new systems of differential equations where, for both A and B configurations, only 4 equations remain coupled in the limit d → 2. For example, in the case of configuration A we can take as new basis , plus the new master defined as where the functions c ij (d; s) are rational functions for the dimension d, the momentum s and the two masses, and are finite as d → 2. This insures, thanks to the overall coefficients (d − 2), that the differential equation for J 5 decouples completely from the other four, as expected.
As far as configuration B is concerned, in order to achieve the complete decoupling of two out of the six equations, we can take as basis (4.60) Using these basis one obtains a new system of differential equations, where the two equations for J B 5 and J B 6 develop an explicit overall factor (d − 2), such that, at every oder in the Laurent expansion, they can be solved trivially by quadrature. Order by order, once the result for the latter is known, one is left with a system of four coupled differential equations for the remaining master integrals. For compactness we prefer not to give here explicitly the systems of differential equations.

The case of three different masses
Generalising even further we can check what happens if three out of the four masses are allowed to take different values. In this case there is of course only one possibility, which we choose to be I 3 (d; n 1 , n 2 , n 3 , n 4 , n 5 , n 6 , n 7 , n 8 , n 9 ) = I 4 (d; n 1 , n 2 , n 3 , n 4 , n 5 , n 6 , n 7 , n 8 , n 9 ) m4=m3 . We start, as always, performing a reduction for generic d. The complexity increases and we find 8 independent MIs We can then consider the system of IBPs for d = 2. It is easy to check that in this case 3 MIs degenerate, and therefore only 5 MIs remain linearly independent. We do not report here the equivalent of relations (4.53) (4.54) and (4.55), since they are considerable more lengthy, but one can easily work out the reduction in d = 2 and find that, for example, I 6 (2; s), I 7 (2; s) and I 8 (2; s) can be written as linear combinations of the I 1 (2; s),...,I 5 (2; s). Using the methods described above, 3 out of the 8 differential equations for this particular mass configuration can be decoupled in the limit d → 2, and one can in principle derive a fifth-order differential equation for any of the MIs, and in particular for the scalar amplitude I 1 (d; s).

The general case of four different masses
Last but not least, we move to considering the most general configuration with four different masses. In this case the complexity increases even further and solving the IBPs in d dimensions brings to a reduction in terms of 11 different MIs The number of independent master integrals is obviously very large and, if all differential equations for the 11 MIs were to be coupled, this would imply an 11-th order differential equation for any of the masters and in particular for the scalar amplitude I 1 (d; s). It is therefore very interesting in this case to know how many MIs can be decoupled using the methods described above. Again it is enough to repeat the reduction to MIs, but fixing this time d = 2, and we immediately find that 5 out of the 11 MIs become linearly dependent and can be expressed in terms of the other 6. Which MIs survive depends of course on the internal algorithm for the solution of the IBPs. In our case we find as independent MIs This implies that, the new basis of 11 d dimensional MIs defined following the recipe given above, fulfills a system of 11 differential equations, 5 of which decouple from the system as d → 2. In this way we expect that a sixth-order differential equation for the scalar amplitude can be derived in d = 2.

Comments and open questions
Before moving on to the conclusions we would like to bring attention to some issues that might have gone unnoticed and which nevertheless leave room to very interesting open questions. In the previous sections we have worked out different applications of the ideas outlined in Section 3. We have seen explicitly that studying the IBPs in d = 2 or d = 4 can provide, for different Feynman graphs, identities useful to decouple the system of differential equations they fulfil. However in this discussion there is a point that we have avoided mentioning on purpose. Let us imagine to have to deal with a Feynman graph with three master integrals I 1 , I 2 and I 3 , which fulfil a system of three coupled differential equations in the limit d → 4, and let us suppose to apply the methods described above in order to try and decouple the system. We can imagine that by solving the IBPs in d = 4 only one relation can be found. In this case we know that we can use it in order to decouple one of the three integrals in the limit d → 4, leaving therefore a system of two coupled equations, equivalent to a second order differential equation. Is this enough to say that there must be no other way to fully decouple all three differential equations? The answer is, in general, of course no. We have discussed already how the decoupling of the differential equations in any even integer number of dimensions is equivalent to the decoupling of the latter in d = 4. In Appendix B we show explicitly how, if one can find a basis that decouples the differential equations in d = 2 n, a corresponding basis can be constructed which decouples them in d = 4. Assuming that this is true, let us then go back to the problem of the three coupled master integrals. Let us imagine that studying the IBPs in d = 2 two linearly independent relations are found among the three masters and that this allows to completely decouple the system in d = 2. In this case we know that a corresponding basis would have to exists in d = 4 as well, and we could find it with the methods described in Appendix B. On the other hand, if we had not found any new relations in d = 2, we could still decide to try in d = 6, or in d = 8, 10, 12 etc... apparently without an end. Who or what tells us when to stop and, therefore, when we can assert without any doubts that not enough relations can be found, in any even number of dimension, to decouple completely the system? This question is extremely interesting and we unfortunately do not have a conclusive answer to it. For sure in all examples considered so far it has always been enough to study the differential equations in d = 2 and d = 4 only in order to find all needed relations to decouple the system. In those cases where not all equations could be decoupled (see for example Sections 4.3, 4.4, 4.6) an attempt to consider different numbers of dimensions would simply produce no new relations at all, suggesting that there is no way to further simplify the problem, at least in this framework. Of course this does not constitute a mathematical proof in any respect. If these considerations have not brought to a definite answer yet, they nevertheless open the possibility for a different perspective in the way a system of differential equations for master integrals should be studied. We usually tend to think that the only physically relevant results are obtained when studying the system in the limit d → 4. A lot of useful information, though, can be extracted studying the system as d → 2 n and, sometimes, the relations found in this way appear to be independent and, in a sense, complementary. Whether these relations are really independent and how to determine the maximum amount of information that can be extracted by studying the IBPs in fixed integer numbers of dimensions remain open questions for now. It seems however reasonable to think that a more global approach, which allows to study the systems of differential equations in general for any even number of dimensions (and not only in the limit d → 4) could bring a much deeper insight in the structure and properties of the latter.

Conclusions
The method of differential equations has proven to be one of the most effective and promising tools for the evaluation of multi-loop and multi-scale Feynman integrals. The usual procedure consists in reducing all Feynman integrals to a basis of master integrals through integration by parts identities, then derive differential equations satisfied by the master integrals and finally try and solve them as Laurent expansion in (d − 4). For many problems of physical interests the coefficients of the Laurent expansion of the master integrals can be expressed in terms of a particular class of special functions called multiple polylogarithms. It has been noted that, whenever this is possible, a basis of master integrals can be found such that their differential equations become triangular in the limit d → 4, allowing a simple integration of the differential equations by quadrature. It was moreover conjectured that in all such cases a canonical basis can be found, turning the integration of the differential equations into a straightforward algebraic problem. If a complete set of boundary conditions is also known, the problem can therefore be considered as completely solved. On the other hand, different cases are known where such a simplification cannot be achieved and a minimum number of differential equations remain coupled. Of course in all these cases also a canonical basis (in the original sense introduced in [28]) cannot be found. Whenever this happens, it becomes of crucial importance to be able to determine the minimum number of master integrals which cannot be decoupled from the system. If two or more equations are coupled, in fact, no general technique exists to find a solution and one must resort to different considerations in order to find a complete set of homogeneous solutions, which can then be used in order to build up the inhomogeneous solution using Euler's method of the variation of the constants (see for example [37]). Of course, the larger the number of coupled equations is, the more difficult it becomes finding a complete set of solutions. Reducing the order of the system of differential equations is therefore essential from a practical point of view in order to be able to successfully tackle the problem. The issue is nevertheless interesting also from a more general point of view. Master integrals satisfying higher order differential equations, in fact, cannot usually be expressed in terms of multiple polylogarithms only and a very intensive theoretical effort has been recently devoted to determine the properties of the new special functions required. The most famous example is the two-loop massive sunrise graph. In this case two differential equations remain coupled and therefore the problem amounts to solving a second-order differential equation. It has been recently shown that the solution of the latter can be expressed in terms of a new generalisation of the multiple polylogarithms, called elliptic polylogarithms. Many questions are nevertheless still to be answered. Are elliptic polylogarithms enough for describing all Feynman integrals whose evaluation can be reduced to a second order differential equation? And what about higher order equations? A first step towards an answer to these questions seems therefore to be in a criterion to determine, given a set of master integrals and the system of differential equations they fulfil, the minimum number of differential equations coupled, and therefore the class of special functions required. In this paper we presented a simple idea which proved to be very useful in this respect. We showed in particular that the study of the IBPs for fixed integer values of the space-time dimensions, d = n, can provide the information required for decoupling the differential equations in the limit d → n. Indeed our criterion is, in principle, a sufficient but not a necessary one, in the sense that we did not prove that if no extra relation can be found among the MIs when d = n, then no decoupling is possible for d → n. The criterion has moreover proven to be extremely effective, inasmuch as it provided, in all cases that we considered, relations useful for decoupling some of the differential equations and therefore substantially simplify the problem at hand. It would indeed be extremely interesting to prove whether this criterion is not only a sufficient but also a necessary criterion, checking, for example, whether the number of independent MIs of the three-loop banana graph in d = 2 (see Section 4.6) can be further reduced by any order means, reducing in this way also the maximum degree of the differential equation satisfied by the scalar amplitude. The criterion is moreover extremely simple to apply, since it can be very easily implemented into any existing public or private IBPs reduction code. In this respect, the possibility of pairing the study of IBPs in fixed numbers of space-time dimensions together with the new concept of a (pseudo-)finite basis of MIs, recently introduced in [54], looks particularly promising. The latter, in fact, could potentially provide a way to automatically determine the highest poles developed by different Feynman integrals for different values of the space-time dimensions.
Acknowledgements I am indebted to Ettore Remiddi, whose continuous advice and support were fundamental for the completion of this work. The idea of looking for identities between master integrals for fixed numbers of the space-time dimensions in order to simplify the systems of differential equations was first developed with him in [39]. The generalisation of those ideas presented in this paper have benefited from many interesting discussions with him. I wish to thank Andreas von Manteuffel for the continuous encouragement to follow the ideas developed in this paper and for having allowed me to use the development version of Reduze 2 prior to its publication. Finally I need to thank Pierpaolo Mastrolia for many discussions and useful input during different stages of the project and Thomas Gehrmann for his continuous support and for his comments on the manuscript.

A Comparison with the Schouten Identities
In this Appendix we would like to show explicitly how the methods described in this paper are equivalent to the Schouten Identities introduced in [39]. We will consider again the two-loop sunrise with one massive and two massless propagators (see Section 4.1 for the definitions of the MIs) and try to derive relation (4.14) using the Schouten Identities. The two-loop sunrise graph depends on three independent momenta, the two loop momenta k, l and the external momentum p. With these three momenta we can build up a d-dimensional Schouten polynomial which becomes zero as d → n with n ∈ N and n ≤ 2. Following [39] we start off by considering the quantity ǫ(k, l, p) = ǫ µνρ k µ l ν p ρ (A.1) defined in d = 3 space-time dimensions. Indeed (A.1) is nothing by the Gram determinant of the three vectors k, l, p. By squaring (A.1) in d = 3 we obtain a polynomial P 2 (d; k, l, p) = (ǫ µνρ k µ l ν p ρ ) 2 = k 2 l 2 p 2 − k 2 (l · p) 2 − p 2 (k · l) 2 − l 2 (k · p) 2 + 2(k · l) (k · p) (l · p) . (A. 2) The polynomial was obtained in d = 3 dimensions, but since it contains only scalar products of the three momenta it can be easily analytically continued to d dimensions and regarded as a d dimensional polynomial. By construction the polynomial is zero as d → 2 the full set of MIs is known as Laurent expansion for d → 2 n, with n ∈ N, by using Tarasov-Lee shift identities [41,42] one can reconstruct their Laurent expansion in any other even number of dimensions, and in particular in d = 4. It is therefore clear that, if by any means one can find a basis of MIs whose differential equations are in a convenient form as d → 2 n (triangular form, canonical form, etc), then there must exists a corresponding basis of MIs whose differential equations look identical under the formal substitution (d − 2 n) → (d − 4). In this Appendix we want to show how this is indeed true and that such a basis can be obtained straightforwardly by a repeated use of Tarasov-Lee identities. Let us start considering a topology with N MIs M j (d; x ij ) where j = 1, ..., N and we made explicit the dependence on the dimensions d and on the invariants of the problem x ij = p i · p j . Let us assume that, similarly to the case of the two-loop massive sunrise, Section 4.3, we are able to find a basis of MIs such that the differential equations take a particularly convenient form in the limit d → 2. In particular, in order to simplify the notation, let us assume that the differential equations are linear in d and can be written as which can rephrased in terms of the new MIs as which is exactly what we were looking for, namely a system formally identical to (B.3) under the replacement (d − 2) → (d − 4). While we showed this for a set of differential equations in a very special form (B.3), the considerations explained above are of course valid for any system of differential equations. Given a set of masters M = (M 1 , ..., M N ) and their system of differential equations in matrix form