On a procedure to derive $\epsilon$-factorised differential equations beyond polylogarithms

In this manuscript, we elaborate on a procedure to derive $\epsilon$-factorised differential equations for multi-scale, multi-loop classes of Feynman integrals that evaluate to special functions beyond multiple polylogarithms. We demonstrate the applicability of our approach to diverse classes of problems, by working out $\epsilon$-factorised differential equations for single- and multi-scale problems of increasing complexity. To start we are reconsidering the well-studied equal-mass two-loop sunrise case, and move then to study other elliptic two-, three- and four-point problems depending on multiple different scales. Finally, we showcase how the same approach allows us to obtain $\epsilon$-factorised differential equations also for Feynman integrals that involve geometries beyond a single elliptic curve.


Introduction
Feynman integrals have played a central role in theoretical physics since the early days of application of perturbative Quantum Field Theory to describe natural phenomena. Together with being indispensable building blocks for the calculation of many relevant observables both for low-and high-energy physics, their explicit analytic evaluation has revealed unexpected connections with modern topics in pure mathematics such as number theory and algebraic geometry. This new point of view has in turn helped to clarify how certain general properties of scattering amplitudes can be made manifest at the level of the Feynman integrals used to represent them. Given the complexity of the calculations required to express scattering amplitudes in terms of special functions, any insights on the types of mathematical structures expected to appear can prove crucial to simplify their computation.
Undoubtedly, a very important role in these developments has been played by the so-called differential equation method [1][2][3], which allows one to demonstrate that Feynman integrals satisfy systems of linear differential equations with rational coefficients. The most straightforward way to derive such differential equations is based on the existence of integration-by-parts identities (IBPs) [4,5] among Feynman integrals. IBPs are linear systems of equations with rational coefficients, which can be solved algorithmically [6] and allow one to express many apparently different Feynman integrals in terms of so-called master integrals. Since these master integrals provide in all respects a basis in the vector space of Feynman integrals for the problem considered, an important question becomes whether different choices of bases can simplify their calculation and highlight the structure of the physical quantities that one tries to compute.
An important step in this direction was achieved a bit more than a decade ago with the introduction of the idea of local integrals [7] as natural building blocks to represent scattering amplitudes in the maximally symmetric N=4 Super Yang-Mills (SYM) theory. A crucial property of these integrals is that they can be expressed as iterated integrations over logarithmic differential forms, with maximum codimension residues normalised to one (or in general to an integer number). These residues are usually referred to as leading singularities and integrals with the property above are said to have unit leading singularities. Interestingly, it was soon realised that under some conditions, these integrals satisfy particularly simple systems of differential equations which are said to be in canonical form [8]. An important (but not the only) property of such canonical differential equations, is that they are in ϵ-factorised form, where ϵ is the dimensional regularisation parameter. Since the whole dependence on ϵ is factorised in front of an ϵ-independent matrix of differential forms, this makes the all-order iterative structure of the Feynman integrals as Laurent series in ϵ manifest. The price to pay to achieve this structure is that the differential forms are no longer guaranteed to be rational functions, but become instead in general transcendental functions. It is easy to see that this is the case since the new forms originate from the solution of the homogeneous system of differential equations close to ϵ = 0, which in turn can be evaluated using generalised unitarity cuts [9][10][11]. A particularly powerful technique to compute these cuts involves the use of the Baikov representation [12][13][14][15].
As long as we limit ourselves to consider integrals that can be expressed as iterated integrations over logarithmic differential forms, a lot is understood about how such a canonical form can be achieved. Geometrically, these integrals can often be related to iterated integrations of rational functions defined on a genus zero surface, i.e. on the Riemann sphere, and in those cases, they can be expressed in terms of a well-understood class of functions, namely multiple polylogarithms (MPLs) [16][17][18][19][20]. 1 Various approaches have been developed to achieve a canonical form in these cases. These include algorithms that aim to transform directly the system of differential equations [22][23][24][25], or that are based on the study of the residues of the corresponding integrands [26]. While none of these techniques is guaranteed to work in all cases, the problem is considered to be at least conceptually under control.
The situation is very different once we leave the space of logarithmic differential forms, increasing the genus and (or) the dimension of the algebraic surfaces involved. This is an important problem, both from a conceptual and a practical point of view, since many examples are known of Feynman integrals which can be naturally defined either on higher-1 Note that this is not always the case, see for example [21]. genus two-dimensional complex surfaces, as elliptic curves, or even on higher-dimensional Calabi-Yau geometries, see for example [11,.
Much progress has been made in the past decade, in particular, for what concerns singlescale problems and the genus one case [11,[48][49][50][51][52]. Still, a lot remains to be understood if either the geometry becomes more complicated or multiple scales are involved. 2 Moreover, also in those cases where results are known, there is not yet a consensus on whether the bases identified fulfil the right criteria to be the natural generalisation of a canonical basis, as originally defined in [7,8]. From this perspective, one of the most pressing questions concerns the possibility to develop an approach that allows one to determine bases of Feynman integrals whose differential equations are in ϵ-factorised-form and that make their iterative structure manifest, for example avoiding linearly dependent differential forms, as much as possible independently of the geometry and of the number of scales involved.
In this paper we try to address this intricate problem, moving from an empirical observation made in the context of Feynman integrals which can be evaluated in terms of elliptic Multiple Polylogarithms (eMPLs) [54][55][56][57]. In fact, it was observed that whenever explicit results can be obtained in terms of a pure version of eMPLs [58], linear combinations of Feynman integrals can be identified that appear to have the right properties of purity and uniform transcendental weight expected from canonical Feynman integrals [59]. 3 This rotation involves a specific separation of the period matrix (which solves the corresponding homogeneous differential equations [9]) into a semi-simple and a unipotent part. This observation has more recently been confirmed in the case of the differential equations of the two-loop sunrise graph [60,61]. The goal of this paper is to show that, if one starts from a basis of master integrals that fulfils some minimal, physically motivated requirements, this construction turns out to be very general and can be used as the starting point towards a procedure to find ϵ-factorised bases of master integrals in a large variety of problems, both in the elliptic case and beyond.
The paper is organised as follows: In section 2 we explain our procedure. Our method is based on five main steps which are introduced and discussed. Subsequently, we provide a pedagogical example of our method in section 3, where we discuss our approach applied to the sunrise family with different mass configurations. We then continue in section 4 with examples of Feynman graph families with a single underlying elliptic curve. We consider one-, two-and three-parameter families of integrals, including a four-point example at two loops. Afterwards, in section 5 we apply our procedure also to Feynman graph families characterised by more complicated geometries than a single elliptic curve. Finally, we provide our conclusions in section 6. We also included several appendices as well as an ancillary file containing explicit results.

Description of the procedure
We consider dimensionally regularised Feynman integral families of the form where ν j ≤ 0 for all j > n. The set of propagators in the denominator {D 1 , D 2 , . . . , D n } is specified by the topology of the Feynman graph under consideration. The {N 1 , N 2 , . . . , N m } are a minimal set of irreducible scalar products in the problem, i.e. scalar products involving loop momenta that cannot be written as a linear combination of the propagators. We set d 0 = 2n for different n ∈ N, depending on the family under consideration. As outlined in the introduction, Feynman integrals that belong to a given family, exhibit the structure of a finite-dimensional vector space, whose basis we refer to as master integrals. We indicate them in vector form as I. Moreover, any choice of basis of master integrals I satisfies a system of linear first-order differential equations with respect to the kinematic variables and the masses the problem depends on. We indicate the set of all variables with z. The system of differential equations is also often referred to as Gauss-Manin (GM) differential system and it takes the general form Our goal is to describe a procedure to find a transformation matrix R(z, ϵ) constructed as a series of subsequent rotations, i.e.
3) that cast the system of differential equations into ϵ-form, namely The crucial property of eq. (2.4) is that the new matrix GM ϵ (z) does not depend on ϵ. If the Gauss-Manin system is in this particular form we call it to be ϵ-factorised. In addition, in the polylogarithmic case one can (conjecturally) always bring the Gauss-Manin matrix GM ϵ (z) in an even more specific form, where all its entries [GM ϵ (z)] ij are given in terms of d log-forms. If this is the case, we say that the system is in canonical form.
An ϵ-factorised form implies that the iterative structure of the solution is manifest in the differential forms appearing in GM ϵ (z) and, if one understands all non-trivial relations among these forms, one can claim to fully control the functional relations among the iterated integrals that stem from them. This statement can be made more precise as follows: If we call the entries of the matrix [GM ϵ (z)] ij = ω ij , then the condition for the iterated integrals to be linearly independent with respect to some subalgebra of functions F, is that there is no exact form η = df with f ∈ F such that ij a ij ω ij = η for a ij ∈ C (2.5) with not all a ij = 0. This condition is then sufficient to guarantee that the corresponding iterated integrals are linearly independent [62]. Clearly, if the forms are all logarithmic, all relations among them descend trivially from corresponding relations among the respective logarithms.
Overview of the procedure Before describing the idea in detail, we would like to stress that our procedure is not fully algorithmic: while we have identified a set of general steps which work very well for all genus one cases we have studied (and also in some cases beyond genus one), the implementation of each step requires a case-by-case analysis and their complexity depends very much on the problem under consideration. We stress that this procedure can be applied to individual sectors of a graph, allowing one to work bottom-up, sector-by-sector, to construct a basis whose Gauss-Manin system is ϵ-factorised. Schematically, we can summarise our procedure in five main steps as follows: 1. Make a reasonable choice of initial basis for the problem considered.
2. Calculate the fundamental matrix of solutions for the homogeneous system of differential equations at ϵ = 0, split it into a unipotent and semi-simple part, and rotate the initial basis with the inverse of the semi-simple part.
3. Adjust factors of ϵ and swap integrals in the basis such that non ϵ-factorised terms only appear below the diagonal of the Gauss-Manin connection matrix. We reach in this way an upper-triangular ϵ-form.
4. Integrate out iteratively the remaining unwanted terms in the homogeneous part of the differential equations.
5. Perform shifts in the inhomogeneous part to factorise ϵ in the whole system of differential equations.
Let us begin by presenting some general considerations which are behind each of the points above.
Step 1: Choice of the initial basis The first step is extremely important but it is also the least understood one. First of all, we should stress the obvious fact that, as long as we assume that a basis of master integrals exists whose differential equations are in ϵ-factorised form, it is always possible, at least in principle, to reach it starting from any other basis by a sequence of transformations involving arbitrarily complicated functions of ϵ and of the kinematics. 4 As a consequence, there is no unique choice of starting basis for our procedure to work. However, as in the purely polylogarithmic case, if we start from a basis that is particularly far from a good basis, our procedure might either become computationally extremely demanding or even fail at any of the subsequent steps. The question becomes therefore, what are general criteria that allow us to say that the basis we are starting from is not too far from a good basis.
To this aim, it is useful to start from the much better understood polylogarithmic case. For a rather large class of problems that can be solved in terms of algebraic functions and Chen iterated integrals over d log-forms, a study of the integrand can provide important information. 5 In fact, as elucidated in [7,8,26], candidates for canonical integrals which fulfil ϵ-factorised differential equations, can be identified by selecting integrands which can be expressed as iterated d log-forms with coefficients equal to numbers. The integrand in this case is said to be in d log-form and to have unit leading singularities. This can be achieved in practice by an analysis of the iterated residues of the corresponding integrands, in any suitable representation (Feynman-Schwinger parameters, Baikov representation, etc). The analysis is usually performed in d = d 0 −2ϵ dimensions, with typically d 0 = 2, 4, 6. As a matter of fact, studying the integrand exactly in d = d 0 provides often enough information to determine a suitable basis for general values of ϵ. 6 Intuitively, this can be understood realising that for d = d 0 − 2ϵ the integrand can be schematically parameterised as with F(x i , z) and G(x i , z) algebraic functions, x i are the integration variables (for example the Baikov parameters) and z the set of kinematical invariants and masses that the integral depends on. Therefore, if the integrand is in d log-form for ϵ = 0, higher order corrections in ϵ will not invalidate this form, as they will naturally add only powers of logarithms. This analysis is extremely powerful and one might wonder how it could be generalised beyond d log-forms. In the following sections, we will provide an example of how this could be achieved in the genus one case, building upon the construction of pure elliptic multiple polylogarithms provided in [58]. Nevertheless, despite being informative, the study of leading singularities has at least two drawbacks. First of all, in a general multi-parameter case and for increasing numbers of loops, it becomes computationally extremely difficult to analyse all iterated residues. In these cases, a simpler analysis restricted to some specific subsets of generalised cuts can provide partial but useful information to define good candidates for a canonical integral. More importantly, depending on the parametrization that one chooses, analysing the residues of the integrand often imposes too restrictive conditions, such that not enough canonical candidates can be identified unless the family of integrals is enlarged. In complicated multi-loop and multi-parameter cases, this is often not practical.
A classical example is provided by integrals with squared propagators. Such integrals are often excluded a priori in the residue analysis since squared propagators typically show up as double poles in the integrand. When dealing with massless Feynman integrals, insisting on the absence of double propagators is often well justified, since the latter typically generate power-like infrared (IR) divergences, which are not expected to appear in gauge 5 Note that this goes beyond the realm of multiple polylogarithms, see for example [21]. 6 One should take extra care if a purely d0-dimensional parametrization of the loop momenta is used, since numerators proportional to Gram determinants would all be identically zero and one would lose possible candidates for canonical integrals. This is avoided using Feynman or Baikov parametrization. theories and which would manifest as non-logarithmic singularities in the Feynman integrals. On the other hand, this analysis is often too superficial, as it does not take into account the fact that, after integration on the contour that defines the Feynman integral, these double poles might rearrange into harmless single poles. This can be seen already in the case of the one-loop bubble integral family: (2.7) It is well known that both for m 2 = 0 or m 2 ̸ = 0, a candidate for a canonical integral in d = 4 − 2ϵ dimensions is provided by the integral with a propagator squared, namely Bub 2,1 (p 2 , m 2 ). The reason why this is a good candidate in both cases can be seen in various ways, the simplest being that the bubble with a squared propagator in d = 4 − 2ϵ dimensions is proportional to the one with linear propagators in d = 2 − 2ϵ. Importantly, in the massless case, one can always trade the bubble with a squared propagator with a threepoint function, since the two are not independent under integration by parts identities. On the other hand, in the massive case, both master integrals are independent and have to be kept. While this example is extremely simple, it points to the fact that admitting bases of master integrals with squared massive propagators is unavoidable to construct candidates for ϵ-factorised differential equations, even in the genus zero case. Moving from these general considerations, we can now summarise our approach to select a good starting point for our procedure. We work sector by sector to construct our basis and all integrals are assumed to live in d = d 0 −2ϵ dimensions, with d 0 = 2n and n ∈ N. We notice here that our procedure can be seen as a generalisation of the one proposed in [63] beyond the polylogarithmic case.
1. First and foremost, we always avoid integrals with power-like UV or IR divergences in d = d 0 .
2. For sectors with one single master integral, we try to select a candidate with unit leading singularities, at least on the maximal cut. In general, more complicated geometries can be hidden in lower cuts. Since we work bottom-up, we assume that an appropriate basis for them has already been worked out.
3. For sectors with more than one master, we first focus on the homogeneous part of the differential equations in strictly d = d 0 . Either by studying the maximal cut of the master integrals in the sector or, equivalently, by analysing the Picard-Fuchs operator associated with the homogeneous system, we determine the geometry associated with that sector (genus zero, genus one, Calabi-Yau etc.). This informs us on the minimal irreducible complexity in the corresponding block of differential equations, i.e., how many master integrals can be decoupled by just rational or algebraic transformations.
4. If the analysis at point 3. reveals that the sector can be completely decoupled, we attempt to choose integrals with unit leading singularities, following the standard approach.
5. If the analysis at point 3. reveals instead that n integrals are coupled even in d = d 0 and if the sector contains m > n integrals in d = d 0 −2ϵ, we first choose the remaining m − n integrals to make this minimally coupled system manifest. This can be often achieved by selecting integrals which are zero in d = d 0 , by introducing numerators build from Gram determinants of loop and external momenta, see for example [64][65][66]. In addition, we also select as many of the extra master integrals as possible, such that their maximal cut can be still localised by taking residues. In the elliptic case, this corresponds to selecting differentials of the third kind.
6. After decoupling, we are left with choosing a basis for the minimally coupled sectors.
In this case, we proceed as follows: First (if necessary, restricting the analysis to the maximal cut) we always choose the first integral to be expressed as a series of d logforms up to at most one integration. For the latter, we impose that it corresponds to the holomorphic differential of the first kind on the geometry considered. In the elliptic case, integrating this form on the two independent cycles provides the two periods of the elliptic curve. Similarly, on more complicated geometries, one obtains one set of independent homogeneous solutions for the corresponding differential equations. One reasoning behind choosing one integral in this form is that integrating over the holomorphic differential form does not introduce any power-like IR divergences. Moreover, as it was worked out explicitly in the elliptic case, this form is one of the integration kernels that we always expect in generalisations of polylogarithms beyond genus zero [58]. Notice that this differential form does not need to show up in the last integration.
7. The last problem is how to choose the remaining integrals in the coupled sectors, since only excluding power-like IR or UV behaviour usually does not completely constrain the basis choice. In the elliptic case, one possibility could be to choose as extra candidates, integrands which can be written as strings of d logs and ∂E 4 ( m n ; x, ⃗ a) /∂x, i.e., kernels that define pure elliptic polylogarithms, see [58]. This can be attempted in simple examples, but the analysis required quickly becomes cumbersome for multiloop, multi-scale problems. Therefore, we choose the extra masters as linear combinations of iterated derivatives of the first integral with respect to the internal masses. This has the advantage of producing differential equations in standard form, without increasing the degree of divergence in the IR. 8. Finally, as a rule of thumb based on experience, we always avoid bases of integrals that generate differential equations whose coefficients have a non-trivial dependence on ϵ in the denominators. 7 Most of the time, imposing the criteria at points 1. to 7. already guarantees that this requirement is satisfied, modulo trivial rescaling in ϵ of the basis elements. If that is not the case, and we still have freedom in choosing our basis, we do that by trial and error avoiding such poles as they increase the complexity of the next steps drastically.
In conclusion, the points above provide us with a guide to constructing a starting basis of master integrals which does not have undesirable properties, among which the most important ones are power-like divergences and non-minimal couplings in the homogeneous blocks of the individual sectors, in the limit d = d 0 . We do not claim that all these points are new elements introduced in this paper but are instead the result of the experience collected by many research groups working on this topic.
Step 2: Rotation by the inverse of the semi-simple part of the period matrix The second step is the crucial one in our procedure. For a choice of basis with no poles in ϵ in the Gauss-Manin connection matrices, we start by computing the fundamental matrix of solutions W at ϵ = 0 for every coupled block in the sector under consideration, discarding any contributions from integrals whose own differential equations do not couple to the block. In the following, we refer to this system also as the maximal cut system, as the maximal cuts of the integrals provide a solution to it [9-11, 14, 15], and to W also as the Wronskian matrix or period matrix. Beyond the polylogarithmic case, W contains new classes of transcendental functions for which we can always derive a representation in terms of locally convergent power series, containing also logarithmic contributions in the parameters. Furthermore, we note that in the construction of W, it is convenient to order the solutions such that the powers of logarithms appearing in their power series expansions increase from left to right in the first row. This is usually referred to as a Frobenius basis.
To proceed, we take inspiration from the patterns observed in [58,59], where pure and uniform transcendental weight elliptic Feynman integrals could be selected by a rotation of the basis of integrals by the so-called semi-simple part of the period matrix. In those references, the observation was based on explicit results obtained by direct integration and therefore was limited to at most two orders in ϵ. Our goal here is to generalise this procedure to all orders in ϵ, working at the level of the differential equations. We then split W into a semi-simple part W ss and unipotent part W u , i.e. (2.8) In general, this splitting is not unique. The only requirements are that the unipotent part W u has to satisfy a unipotent differential equation [67], i.e. of the type where the matrices U i (z) are nilpotent matrices, while the semi-simple matrix should only be invertible. For the purpose of our procedure, we perform the splitting in such a way that the semi-simple part W ss has lower triangular form, while the unipotent piece W u has upper triangular form, with its diagonal containing only constant entries, normalised to one for convenience. We expect such a splitting to always be possible for a choice of initial basis following the guidelines outlined in the previous paragraph. We then rotate the basis in the top sector with the inverse of W ss , discarding the contributions of the solutions contained in the unipotent part.
Step 3: Transformation to upper triangular ϵ-form The third step of the procedure involves typically nothing more than just adjusting some factors of ϵ, as well as potentially swapping the position of some of the integrals in the basis, such that terms that are not ϵ-factorised appear only below the diagonal of the new Gauss-Manin connection matrix GM. Further, all of these terms are either constant or go with inverse powers of ϵ. Notice that at this point, poles in ϵ typically appear.
Step 4: Clean up of the homogeneous differential equations As the next step, we integrate out these remaining terms systematically in each homogeneous block. This is achieved by shifting the master integrals which generate non ϵfactorised homogeneous equations, by terms proportional to the other master integrals in the same sector. Frequently, this amounts to removing total derivatives involving the objects already introduced by the rotation by the inverse of W ss . However, for more complex problems it happens that not all non-ϵ-factorised terms can be removed in this way. In these cases, it becomes necessary to introduce new functions, which are defined as (iterated) integrals of products of rational or algebraic functions and the transcendental functions introduced by W ss . This leads to new types of integration kernels in the final ϵ-factorised differential equations. It is desirable that these new objects are linearly independent under integration by parts identities as discussed around eq. (2.5). We will comment on this below when looking at explicit cases.
Step 5: Clean up of the inhomogeneous differential equations After having completed the ϵ-factorisation of all homogeneous blocks, the factorisation of the whole system can be achieved by performing suitable shifts in the definitions of the master integrals of a given sector, by integrals in lower sectors. This is in principle a straightforward procedure. In practice, however, it can become tricky as also the functional dependence on ϵ of the coefficient functions in these shifts might be complicated. Moreover, if a given subsector exhibits singularities which were not present in the homogeneous equations for the sector considered, the introduction of new functions might be required. These are likewise (iterated) integrals of functions already present in the differential equations.
3 Application to the two-loop sunrise family Let us see explicitly how our procedure works in the case of the two-loop sunrise graph [28][29][30][31][68][69][70][71][72][73][74][75] shown in figure 1. We work as it is customary in d = 2−2ϵ space-time dimensions and, as it is well known, the result in d = 4 − 2ϵ can be obtained by dimensional shift [76]. We take the associated integral family in the general case of different masses to be defined by the following three propagators and two irreducible numerators: We assume m 1 , m 2 , m 3 > 0 and, in the remainder of this section, we set p 2 = s. As the first step, we analyse the integrand of the so-called corner integral I 1,1,1,0,0 to determine the complexity of the problem. We use a loop-by-loop Baikov representation [12][13][14][15] in two space-time dimensions. Working out the change of variables we find From this representation, it is easy to see that the maximal cut is proportional to a one-fold integral over the reciprocal of the square-root of a quartic polynomial in the integration variable z 4 As long as none of the internal masses vanishes, all four roots are distinct and the geometry is that of an elliptic curve. If at least one of the internal masses is zero, two of the roots degenerate and the elliptic curve reduces to a genus zero surface. In the following, we consider three cases explicitly: 1. Two equal non-vanishing masses 0 ̸ = m 2 2 = m 2 3 =: m 2 and one vanishing mass m 2 1 = 0. This case is polylogarithmic and we can see how our approach works in the case of genus zero geometries.

Three equal non-vanishing masses
This is the simplest genuine elliptic case.
3. Three non-vanishing masses of which two are equal, 0 ̸ = m 2 1 ̸ = m 2 2 = m 2 3 ̸ = 0. This case is multivariate and of similar complexity as the generic case of three different internal masses but with considerably more compact mathematical expressions.

The sunrise with two equal non-vanishing masses and a vanishing mass
For m 2 2 = m 2 3 =: m 2 and m 2 1 = 0, the sunrise integral family admits three master integrals, a tadpole and two in the top sector. In this case, the underlying geometry is a Riemann sphere, so if we were to follow our recipe above, we should attempt to identify integrals of unit leading singularity in d = 2 space-time dimensions. It is immediate to do this analysis directly or using for example DlogBasis [26]. Three candidates are quickly found to be associated to the threshold for the production of two massive and one massless particle.
We chose it such that it is manifestly real above the threshold s ≥ 4m 2 . The factor ϵ 2 is purely conventional such that the expansions of all integrals start at order ϵ 0 . It is easily checked that these master integrals satisfy indeed differential equations in canonical form. Instead of following this standard approach, let us pretend that we were unable to identify all these candidates and their normalisation factors in the above analysis and use this example to illustrate the steps involved in our procedure. As we will see, they yield essentially the same canonical basis. We start from the basis of master integrals I = {I 1 , I 2 , I 3 } given by They were chosen based only on the following knowledge of their properties: 1. The tadpole integral I 1 is trivial and known to evaluate to a pure function of uniform transcendental weight in d = d − 2ϵ.
2. The first integral in the coupled sector, I 2 , has a representation in terms of iterated d log forms in d = 2 space-time dimensions.
3. The second master, I 3 , is proportional to the derivative of the first one with respect to the internal mass squared. In practice, this means that we put a dot on one of its two massive propagators. This guarantees that no spurious power-like IR divergences are introduced.
The system of non-ϵ-factorised differential equations for our starting basis I reads dI = GM m 2 dm 2 + GM s ds I , For simplicity, we work with the Gauss-Manin connection matrix with respect to m 2 . The one with respect to s follows at each step from a scaling relation implied by dimensional analysis and Euler's theorem on homogeneous functions 8 . Concretely, we have in (3.10) The next step in our procedure is to analyse the maximal cut system, which corresponds to the following set of differential equations ∂ ∂m 2 (3.12) Its fundamental matrix of solutions W can be chosen to take the following form 9 where r(s, m 2 ) was defined in eq. (3.8). In principle, it is possible to rotate the basis of the top sector with the inverse of W to solve the system from eq. (3.12) at ϵ = 0. However, this step does not lead to a canonical basis, not even to a factorisation of ϵ. Moreover, the coefficient functions in the differential equations will be of mixed transcendental weight. Instead, we follow Step 2 of our procedure and perform the splitting of W into a lower-triangular semi-simple part W ss and an upper-triangular unipotent part W u . This works because our choice for the first integral in the top sector I 2 already has uniform transcendental weight as we anticipated from the properties of its integrand. Explicitly, we find such that W u satisfies the unipotent differential equation Alternatively, one might also reduce this to a one-variable problem by working with the ratio of the two scales, s/m 2 . 9 If we had started with a better choice of initial master integrals I, i.e. I1,1,1,−1,0 instead of I1,1,2,0,0, the Wronskian matrix (3.13) would have been simpler. Here, we choose a non-optimal basis intentionally to illustrate more features of the procedure.
By splitting the Wronskian matrix W in this way, the logarithm appears exclusively in the unipotent part, while the semi-simple part has uniform transcendental weight zero. Rotating now the basis in the topsector only with the inverse of W ss , we arrive at a new basis which satisfies the differential equation (3.17) The remaining steps to factor out ϵ in this problem are now straightforward. First, we rescale the last integral in the basis by a factor 1/ϵ with respect to the other two. An additional factor ϵ 2 is added to all integrals for the same conventional reasons as for the basis in eq. (3.7). This leaves us with just a single term that is not yet proportional to ϵ. It is, however, a total derivative of an algebraic function already appearing in the problem and it is therefore easily integrated out. These manipulations can be summarised in the rotations which allows us to define a basis that satisfies canonical differential equations As advertised, the relation to the canonical master integrals found by DlogBasis is a simple constant rotation that is given by We would like to stress the following points concerning the application of our procedure in the polylogarithmic case: • The critical step was the splitting of the Wronskian in the semi-simple and unipotent part. Discarding the latter, allowed us to remove the logarithm appearing in the Wronksian matrix, which is the function that would have ruined the transcendental weight properties and the iterative structure of the differential equations.
• In general, we do not expect our procedure to outperform other approaches to find canonical bases in the polylogarithmic case. In particular, for an increase in the number of master integrals, the steps involved also increase in complexity. However, if analysing the leading singularities does not provide enough candidates, our method could be used complementary.

The sunrise with three equal non-vanishing masses
We now take one step up in complexity to the simplest elliptic scenario and set all three internal masses of the sunrise to the same value, i.e. m 2 1 = m 2 2 = m 2 3 =: m 2 . The number and sector distribution of master integrals does not change compared to the previous case.
First, we need to choose a basis of integrals to start from. Taking the limit for equal masses in eq. (3.3), it is easy to study all iterated residues of the integrand corresponding to I 1,1,1,0,0 . The complexity of this analysis depends in general on the order in which the various residues are taken. For this particular problem, it is convenient to take residues in the natural order z 1 to z 4 . For the first three integration variables, one easily sees that the integrand has only simple poles in the corresponding variable, such that as expected the result can be written schematically as follows: where the f i are algebraic functions, at most quadratic in the integration variables z 1 , z 2 and z 3 . If the last integration in z 4 were missing, this integral would be a series of iterated integrations over d log forms and would therefore have the right characteristics to be a candidate canonical integral in the language introduced in [7]. Interestingly, as already hinted at in the previous section, eq. (3.21) shows that the last integration involves only the holomorphic differential form of the first kind on the elliptic curve, which corresponds exactly to one of the integration kernels that define elliptic multiple polylogarithms, as it is manifest in the form originally proposed in [58]. Elliptic multiple polylogarithms can in fact be defined as iterated integrations of rational functions on an elliptic curve as follows: with n i ∈ Z, c i ∈ C and ⃗ a = {a 1 , ..., a 4 } is the vector of roots defined in eq. (3.6). The first integration kernel for n = 0 reads explicitly where y = P 4 (z 4 ) defines the elliptic curve, ϖ 0 is the holomorphic period on the curve and c 4 is an algebraic function which is added for convenience. With these definitions, and assuming for the sake of the argument that the branching points of the elliptic curve ⃗ a are constant, we can schematically write By generalising the original definition of local integrals provided in [7], we can conjecture that an integral in this form is a good candidate for a "canonical" integral defined on an elliptic curve. We take, therefore, as starting basis I = {I 1 , I 2 , I 3 } defined as where again I 3 has been chosen to be proportional to the derivative with respect to the internal mass m 2 of I 2 . This basis satisfies the system of differential equations where, as before, the Gauss-Manin connection matrix with respect to s follows from a scaling relation. The tadpole integral I 1 is the same as before and remains a suitable candidate for a canonical basis, so we focus our efforts on the top sector. Its homogeneous system of differential equations at ϵ = 0 reads ∂ ∂m 2 A basis for the solution space for I 2 can be constructed explicitly by computing the integral from eq. (3.4) over two independent contours dubbed cycles. To do so, we consider the kinematic region above the threshold for the production of three massive particles, s > 9m 2 , such that the branch points of the elliptic curve (3.6) are ordered according to Two independent cycles can then be chosen as the counter-clockwise contour around the branch cut from a 1 to a 2 and a similar contour encircling the branch points a 2 and a 3 . The resulting two elliptic functions with a convenient normalisation are given by where we have and the complete elliptic integral of the first kind K(k 2 ) is defined by Notice that the first solution is holomorphic in a neighbourhood of m 2 /s = 0, while the expansion of the second one contains log (m 2 /s). Given these, we can now construct the Wronskian matrix to take the following form where we suppressed the dependence of the periods on s, m 2 for ease of typing. We proceed by splitting W into its semi-simple and unipotent part. As outlined before, this splitting is done by demanding the unipotent part to take upper-triangular form with its diagonal normalised to one, while the semi-simple part should be lower-triangular. This yields for the unipotent part where we introduced the τ -parameter defined as τ = 1 2πi (3.33) To write down the semi-simple part, we make use of the Griffith's transversality conditions [77][78][79][80], which in the elliptic case correspond to the well-known Legendre relation for elliptic integrals, The left-hand side of (3.34) is proportional to the determinant of W. Consequently, the function on the right-hand side can in practice also be computed conveniently from Abel's identity as the solution of a first-order differential equation with rational coefficients. The constant of integration depends on the normalisation of ϖ 0 and ϖ 1 , which is, however, irrelevant to our goal. With this formula at hand, the semi-simple part of W can be written as (3.35) Once we reach this point, the remaining steps work as in the polylogarithmic case considered in the previous subsection. After rotating the top sector with the inverse of the semi-simple part, we perform a simple rescaling and integrate out a total derivative to arrive at ϵ-factorised differential equations for the equal-mass sunrise graph: It can be verified that the basis J is, up to a simple constant rotation, the same basis as derived with other methods in [48,81]. The rotation T in (3.36) is constructed from functions with rational dependence on s, m 2 , ϖ 0 and ∂ m 2 ϖ 0 . Notice, however, that the final ϵ-factorised Gauss-Manin connection matrix in (3.36) only contains ϖ 0 and not its derivative. Moreover, all of its entries have at most simple poles in all singular limits, because ϖ 0 ∼ 1/m 2 as m 2 → ∞. If desired, it is also possible to change variables to the canonical variable τ defined in (3.33). Then the entries of GM ϵ m 2 will be given by modular forms of the corresponding monodromy group of the sunrise Γ 1 (6), as it was first observed in [75,82]. On a side note, we would like to mention that upon replacing ϖ 0 by any Qlinear combination of ϖ 0 and ϖ 1 in T would lead to the identical ϵ-factorised Gauss-Manin connection matrix ϵ GM ϵ m 2 with ϖ 0 replaced by the same linear combination. Finally, we could have started from a different basis, where the second master integral in the top sector could be chosen by embedding the sunrise in the kite integral family, and selecting as independent a reducible integral with an extra massless propagator. Starting from this different basis, one can prove that our procedure produces an ϵ-factorised basis that is identical to the one obtained here, modulo a rotation by a constant numerical matrix.

3.3
The sunrise with three non-vanishing masses of which two are equal As a third introductory example, we consider now the fully massive case, where two of the three internal masses are equal but different from the third mass. This increases the number of scales in the problem by one compared to before and gives rise to a second independent tadpole master integral, as well as to a third master integral in the top sector. As initial basis we take the master integrals I = {I 1 , . . . , I 5 } given by The explicit expressions for the Gauss-Manin connections exceed what can be printed in a readable format, but they can be found in an ancillary file.
Let us now comment on how we chose the starting basis in eq. (3.37). The tadpoles are trivial, while we know from the previous example how to choose the first and second integral in the sunrise top sector. From the maximal cut in eq. (3.4), we know that the geometry remains elliptic in even numbers of dimensions. Equivalently, the Picard-Fuchs operator associated to the sunrise sector factorises at ϵ = 0 into a second-order and a firstorder irreducible operator. We expect, therefore, that there should be a minimally coupled system of two master integrals in d = 2, and not one of three master integrals.
To make this decoupling manifest, the integral I 5 was chosen by analysing the integrand associated to the maximal cut in a loop-by-loop Baikov representation as in (3.2) and realising that when the masses are different, one can choose an independent candidate which possesses a non-vanishing, constant residue at z 4 = ∞. Explicitly, we have with P 4 (y) = y 4 P 4 (1/y) . Similarly, one could decouple the third integral using a Gram determinant as it was shown in [64]. The maximal cut system of the (2 × 2)-block formed by the master integrals I 3 and I 4 encodes the elliptic part of the sector for which we want to compute the semi-simple part of the Wronskian matrix. This works exactly as in the equal-mass case. Nevertheless, we will see something new happening in this case. Let us quickly go through the steps involved. Explicitly, the maximal cut system reads where In the kinematic region above the threshold for the production of three massive particles, s > (m 1 + 2m 2 ) 2 , the elliptic functions spanning the solution space can be written as (3.43) We will again omit the explicit dependence of the periods on the kinematic variables whenever convenient. The resulting Wronskian matrix takes the form The unipotent piece is analogous to the equal-mass case but we will not need its explicit form in the following. What we do need is the semi-simple part which, after employing the Legendre relation, reads With this ingredient, we can now construct the following rotation to apply to our initial basis I to reach upper-triangular ϵ-form: Notice that this step includes swapping the fourth and fifth integral in the basis after the rotation coming from the semi-simple part of the Wronskian matrix. Finally, we want to integrate out the remaining non-ϵ-factorised terms. Here, however, one immediately realises that, in contrast to the equal-mass case, they are not just total derivatives of rational functions in s, m 2 1 , m 2 2 , ϖ 0 and ∂ m 2 2 ϖ 0 . To integrate them out, we need to introduce a new function, which we call G := G(s, m 2 1 , m 2 2 ). The new function G can be defined by its partial derivatives (3.47) Using the differential equations for ϖ 0 , it is easy to verify the Schwarz integrability conditions The last line in eq. The lower boundary was set to zero for convenience. We stress, however, that the particular choice of integration constant won't impact the factorisation of the differential equations. Similarly to (3.48), by using the differential equations for ϖ 0 under the integral sign, one can show explicitly that (3.49) also satisfies the other two relations in (3.47). By using the explicit representation for ϖ 0 from eq. Π a, k 2 + 4m 2 2 ϖ 0 + 1 . We recall here that the elliptic integral of the third kind Π(n, k 2 ) is defined as That G admits a representation in terms of an elliptic integral of the third kind is not at all surprising. In fact, if we solved the homogeneous system at ϵ = 0 for the whole 3 × 3 sunrise sector, the corresponding Wronskian matrix would contain terms related to integrating the maximal cut of I 5 (c.f. eq. (3.39)) over the cycles of the elliptic curve. Its integrand corresponds to the differential of the third kind on an elliptic curve and its integration yields elliptic integrals of the third kind, as can be verified readily by an explicit computation. We expect therefore functions of this type to appear in the Laurent expansion of the sunrise master integrals to arbitrary orders in ϵ. In particular, one can verify that when integrated over the cycle of the elliptic curve for which MaxCut(I 1,1,1,0,0 ) ∝ ϖ 0 . After introducing G, the factorisation of ϵ is achieved as follows: (3.54) The terms in eq. (3.54) that are independent of G are required to integrate out all terms containing ∂ m 2 2 ϖ 0 in the Gauss-Manin connection matrix for the variable m 2 2 . After this step, there are remaining terms below its diagonal, that are not yet proportional to ϵ, but that depend on ϖ 0 . Hence, we need to introduce a function whose derivative is proportional to ϖ 0 , from which G emerges naturally. The ϵ-independent Gauss-Manin connection matrices GM ϵ z for z = s, m 2 1 , m 2 2 read explicitly as follows: where we have The P i (s, m 2 1 , m 2 2 ) for i = 1, 2, 3, are polynomials whose explicit expressions read We would like to point out that, if we rescale the last integral in the basis by a factor −2/3, the entries of GM s,ϵ z become persymmetric. Notice that, as in the previous case, none of the kernels in the Gauss-Manin (3.54) contains derivatives of the period ϖ 0 . While this allows us to exclude obvious integration by parts identities among the integration kernels, non-trivial relations cannot be ruled out and should be investigated explicitly case by case.
It is interesting to study how the basis J yielded by our procedure behaves in the limits m 1 → 0 and m 1 → m 2 , i.e. how it is related to the ϵ-factorised bases from the previous two subsections. In both cases, the integrals I 1 and I 5 in our starting basis can be written as linear combinations of I 2 , I 3 and I 4 . Moreover, for s > 9m 2 > 0 one can prove the following functional relations The objects r(·, ·) and ϖ 0 (·, ·) on the right-hand-sides have been defined in eqs. (3.8) and (3.29). The easiest way to verify these relations is from the explicit representations for ϖ 0 and G in terms of elliptic integrals in eq. (3.42) and eq. (3.50), respectively. Using the above relations, it can be checked that all the integrals in the basis J reduce to simple Q-linear combinations of the ϵ-factorised bases we found in the previous two subsections, c.f. eq. (3.19) and eq. (3.36).
Before we conclude this section let us comment on the case of three different internal masses. This setup gives rise to a third independent tadpole, but more importantly, a fourth master in the top sector, for which a convenient choice of initial basis integral is given by I 1,1,1,0,−1 . The rotation to an ϵ-factorised basis is constructed in an analogous fashion as in eqs. (3.46) and (3.54). However, as expected, one now has to introduce two new functions in the last step to integrate out all undesired terms, and which correspond to two differentials of the third kind. These can again be represented as integrals over the corresponding holomorphic period but with additional rational functions in the integrand. The explicit expressions are provided in an ancillary file to this manuscript.

Applications to Feynman integrals with a single elliptic curve
In the previous section, we have worked out explicitly various mass configurations for the two-loop sunrise graph, showing how our procedure allows us to obtain ϵ-factorised systems of differential equations in almost algorithmic steps. For the multi-scale cases, this came at the cost of introducing extra objects, defined as integrals over the holomorphic period of the elliptic curve and rational functions. We have seen that these can be related to differential forms of the third kind. In this section, we will present several examples of single-and multivariate Feynman integral families related to a single elliptic curve, where our approach can be successfully applied. We will see that the patterns we observed can be extended also to three-and four-point functions.

A single-scale elliptic three-point function
As the first example beyond the sunrise topology, we consider the triangle graph shown in figure 2, for which we use the following integral family [59,83] 10 (4.1) We take p 2 1 = p 2 2 = 0 and s = (p 1 +p 2 ) 2 such that the integrals, suitably normalised, depend on the single variable At variance to the sunrise, we consider this two-loop triangle in d = 4 − 2ϵ dimensions and we also set m = 1 for simplicity.
Upon IBP reduction, we find that the problem has 11 master integrals. Following our general prescription, a suitable starting basis I = {I 1 , . . . , I 11 } is given by the integrals  The master integrals {I 1 , . . . , I 9 } are polylogarithmic and their differential equations can be brought into canonical form by standard methods. Explicitly, by studying the leading singularities of the integrals above, one can find the rotation by which a complete ϵ-factorisation of the subsectors can be obtained [83]. Now we can focus on the last two master integrals. As for the sunrise, integral I 10 is chosen by studying the leading singularities of the top sector integrals and selecting an integrand that can be put in the schematic form of eq. (3.24). The candidate for I 11 is instead selected by adding a dot on a massive propagator, which does not worsen its IR behaviour. Moreover, it is convenient to normalise the polynomial defining the underlying elliptic curve (c.f. (3.23)) such that its highest monomial has coefficient one. This implies a normalisation of the master integrals I 10 and I 11 by −1/z. Then we consider their homogeneous differential equations d dz . (4.5) This system is of elliptic type and has as fundamental solutions (4.6) From this we can consider the Wronskian matrix satisfying (4.5). We write its semi-simple part as and as before we rotate our basis of masters by the inverse of W ss followed by a suitable ϵ-rescaling of the integrals. After that, one immediately sees that the term proportional to ϖ ′ 0 (z) and the non-ϵ-factorised term can be removed by shifting the integrals by a total derivative. The full rotation can then be obtained as (4.9) With this, the ϵ-factorised GM equations are given by where the Gauss-Manin connection GM ϵ is given in appendix A eq. (A.1). We see here that in this single-scale problem we do not need additional new functions. In fact, only the holomorphic solution ϖ 0 is needed, while its derivative ϖ ′ 0 does not appear in GM ϵ . The absence of the derivative of the period is an indication that all differential forms are independent under integration by parts identities. Nevertheless, we stress that we have not proven this last statement formally and we cannot exclude that non-trivial relations exist. The next family of Feynman integrals we consider is given by another triangle graph (see figure 3) but with a different mass configuration of the propagators and numerator

A second single-scale elliptic three-point function
(4.11) For the external momenta, we again consider p 2 1 = p 2 2 = 0 and work with the dimensionless variable As before, we study these integrals in d = 4 − 2ϵ and set for simplicity m = 1. We recall here that this class of integrals contributes to the so-called electroweak form factor and was first studied in [32] and later on in [59]. 11 Integration by part identities show that there are 18 independent master integrals I = {I 1 , . . . , I 18 } for which we choose the following starting basis I 1 = I 0,0,0,0,2,2,0 , I 2 = I 0,0,2,1,0,2,0 , I 3 = I 2,0,1,0,0,2,0 , I 4 = I 0,2,1,0,2,0,0 , (4.13) The first fifteen integrals {I 1 , . . . , I 15 } are of polylogarithmic type and a rotation T sub that puts them in canonical form is readily found by standard methods. The explicit form of this transformation is given in appendix A, see eq. (A.2). Let us focus instead on the elliptic top sector, which this time contains three master integrals I 16 , I 17 , I 18 . A leading singularity analysis for these integrals also suggests to normalise them with a factor of z, such that the quartic polynomial defining the underlying elliptic curve has the highest monomial normalised to one. This also has the effect of decoupling completely the integral I 18 in d = 4 from the other two, since the numerator that defines it, generates a residue at infinity exactly equal to 1/z. With this, the maximal cuts of the normalised integralsĨ 16 ,Ĩ 17 satisfy an elliptic homogeneous set of differential equations d dz (4.14) The two fundamental solutions are given by 12 ϖ 0 (z) = 16 π K 16 The semi-simple part of the corresponding Wronskian is given by (4.16) 12 We have neglected terms proportional to log(2) in the expansion of the second elliptic function which is a result of the variable z requiring a rescaling to yield a non-integer coefficient expansion of the first elliptic function. As it was observed in [32] one can also make a change of variables and a subsequent rotation to relate the two elliptic functions (4.15) of this triangle family to the elliptic functions of the equal-mass sunrise graph in two dimensions.
From this we can construct, by the same reasoning as before, the full rotation matrix (4.17) The final ϵ-factorised GM differential equation is then given by where the GM connection GM ϵ can be found in appendix A eq. (A.3). We stress that once again, ϖ ′ 0 does not appear in GM ϵ , an important prerequisite to conjecture that all differential forms in eq. (4.18) are linearly independent on the field of algebraic functions.
As for the equal-mass sunrise and the first triangle family studied in the previous subsection, at every step of the rotation to an ϵ-factorised basis we need only ϖ 0 and ϖ ′ 0 and no extra function has to be introduced. This might be surprising since, contrary to the two previous cases, there is a third master integral in the top sector, whose integrand on the maximal cut at ϵ = 0 corresponds to a differential of the third kind on an elliptic curve. However, upon explicit integration over the cycles of the elliptic curve, the resulting elliptic integral of the third kind might degenerate to linear combinations of ϖ 0 , ϖ 1 and their derivatives, which would explain why no new function is required to achieve the factorisation of ϵ. For the region 0 < z < 1, we verified that upon integration over the cycle of the elliptic curve for which MaxCut(Ĩ 16 ) ∝ ϖ 0 , one finds indeed such a relation: where c is just a constant.

Two-parameter triangle
We have seen that in the different mass sunrise graph, which was also the only multiparameter problem analysed so far, new objects had to be introduced to achieve a complete ϵ-factorisation of the differential equations. In this and the next subsection, we study other multi-parameter problems, related to more complicated processes, to see if the same pattern can be observed. As a two-parameter family, we consider again the non-planar triangle from figure 2, but this time with one more external leg off-shell, i.e. p 2 1 = 0 and p 2 2 = M 2 . The problem depends now on two dimensionless parameters, which we choose to be z 1 = s/m 2 and z 2 = M 2 /m 2 .
As before, we hide the dependence on the variables z 1 , z 2 in the periods and abbreviate ϖ 0 = ϖ 0 (z 1 , z 2 ) whenever convenient. Finally, we removed the terms which are not proportional to ϵ by simple rotations.
To achieve a complete ϵ-factorised form for the full Gauss-Manin system, we have to perform additional rotations in the subsectors contributing to the inhomogeneous part of the differential equation for the top sector. Compared to the single-scale problem analysed before, in this case, not all non-ϵ-factorised terms can be removed by simple rotations that do not involve any new function. There are in fact two terms which can only be removed by introducing a new object, say G := G(z 1 , z 2 ), defined by the two partial derivatives One can easily see that the first equation is solved by The second equation in (4.24) is then satisfied by using the differential equations for ϖ 0 . We stress here that, at variance with the two-parameter sunrise case, where a new function G was required for the ϵ-factorisation of the (3 × 3)-homogeneous system, in this case we do not need a new function to obtain an ϵ-factorised form for the homogeneous system associated with the top sector. Instead, the new function originates from the subtopologies, which can be understood as follows: First, the additional scale that this problem depends on does not increase the number of master integrals in the top sector. This can be interpreted as the fact that, contrary to the sunrise, there is no extra residue in the integrand of the maximal cut, which is linearly independent under integration by parts. This implies that there is no contribution from an elliptic integral of the third kind to the homogeneous differential equations. Second, there is instead a different singularity structure in the subsectors, namely a new pole compared to the ones found in the homogeneous differential equations for the top sector. Integrating over this pole gives a contribution formally similar to an extra residue in the homogeneous equation. The full rotation matrix and the final form of the GM matrices can be found in an ancillary file.

Three-parameter double box
As a final example, we consider the double box graph depicted in figure 4. The associated integral family is defined by the following propagators and irreducible numerators: The external momenta satisfy p 2 1 = p 2 2 = p 2 3 = 0 and p 2 4 = M 2 . With Mandelstam variables defined by s = (p 1 + p 2 ) 2 and t = (p 1 + p 3 ) 2 , the family depends on a total of four dimensionful or, by dimensional analysis, three dimensionless variables. Further, we set d = 4 − 2ϵ for this example.
This graph contributes to the planar corrections to the production of a Higgs boson and a jet through a loop of massive quarks and was studied at length in the literature, see [85,86], where a complete calculation could only be achieved numerically. 13 In the literature, a total of 73 master integrals had been identified. A new reduction performed with Kira 2.0 [87] reveals an additional relation, originating from a higher sector, which reduces this number to 72. The relation reads (4.27) In [85], 64 independent candidate integrals for a canonical basis for all but one subsector had been identified (in the notation of [85], f A 19 is no longer independent and we remove it from our basis). Moreover, four suitable candidates on the maximal cut of the top topology have been derived. The remaining four integrals lie in an elliptic subsector depicted in figure 5, which is obtained from the top topology by pinching line number three. Here, we wish to complement these results by applying our method to find suitable basis integrals for the elliptic sector, such that a complete basis satisfying ϵ-factorised differential equations is at hand. In the following, we will give an overview of the relevant details of each step. As a starting basis to apply our method, we use a slight variation of the integrals 13 Note that compared to [85], our propagator definitions (4.26) differ by a sign.
• The maximal cut of I 1,1,0,2,1,1,1,0,0 is proportional to ϵ which means that it vanishes in strictly d = 4 space-time dimensions. Therefore, its differential equation at ϵ = 0 should not couple to I 1,1,0,1,1,1,1,0,0 , I 2,1,0,1,1,1,1,0,0 and I 1,1,0,1,1,1,1,0,−1 , as this would imply a non-trivial relation among their maximal cuts. However, this also means that rescaling it by a factor 1/ϵ with respect to the other three master integrals makes it a potential candidate for our starting basis and we can verify from the explicit form of the differential equations that it is indeed suitable. Further, the maximal cut of I 1,1,0,2,1,1,1,0,0 possesses a non-vanishing residue, in contrast to the other three master integrals. From this one can read out the normalisation factor (4.29) which has to be added to the transformation to achieve full ϵ-factorisation.
• The maximal cut of the first integral I 65 at ϵ = 0 can be written in a form analogous to eq. (3.24), but where the elliptic polylogarithmic kernel does not appear in the last integration but in the next to last integration. As discussed above, conjecturally this is still a good candidate for one of the master integrals.
• Finally, the fourth integral I 66 was chosen with a dot on a massive propagator such that it is independent of the other three. This excluded putting the dot on the seventh propagator.
All subsequent steps to achieve the factorisation of ϵ proceed as in previous cases, without any additional complications. The elliptic function introduced by the rotation with the inverse of the semi-simple part of the Wronskian solution matrix can be chosen as (4.30) Following the discussion in the previous section, in this case, as expected from the number of independent master integrals in the top sector, the ϵ-factorisation of the homogeneous system requires the introduction of two new functions G 1 , G 2 . In addition, two new functions G 3 , G 4 are required to factorise ϵ in the whole system, including the subtopologies. These can be written as (4.32) The square-root in the denominator of r 2 (s, t, M 2 , m 2 ) stems from the normalisation factor (4.29) required for the integral I 67 in our starting basis. The expressions for r 3 (s, t, M 2 , m 2 ) and r 4 (s, t, M 2 , m 2 ) contain additional square-roots, generated by other sectors. Thereby, we would like to point out that the function G 4 is required to factorise ϵ in the contributions from the elliptic sector to the differential equations of the top sector.
Also in this case, as for all previous examples, the final ϵ-factorised Gauss-Manin connection matrix does not contain any explicit dependence on derivatives of ϖ 0 . The explicit rotations and expressions for the Gauss-Manin connection matrices are included in an ancillary file.

Cases beyond a single elliptic curve
In the previous sections, we have provided examples of how our procedure can be applied to multi-scale elliptic problems to obtain fully ϵ-factorised systems of differential equations. Here, we want to test the applicability of our procedure beyond the elliptic case, considering families of Feynman graphs with different underlying geometries. We start with a family with two elliptic curves and then consider a three-loop example characterised by a K3 surface.

The three-loop ice cone
Our first example of a Feynman graph family with underlying geometry beyond a single elliptic curve is the three-loop ice cone family (see figure 6). Figure 6. The three-loop ice cone graph.
As it was argued in [66], by studying the maximal cut of the three-loop ice cone in d = 2 one finds two different elliptic curves, which can be both related to the elliptic curve of the two-loop sunrise graph. We follow the conventions from [66] and define the propagators and irreducible scalar products of the ice cone family as Again following [66], we take as the set of starting master integrals 14 I 1 = I 1,1,1,0,0,0,0 , I 2 = I 1,1,1,1,0,0,0 , I 3 = I 2,2,0,1,1,0,0 , I 4 = I 1,1,1,1,1,0,0 , I 5 = I 1,1,1,1,1,−1,0 , I 6 = I 2,1,1,1,1,0,0 , I 7 = I 2,1,1,1,1,−1,0 , To simplify the analysis of the leading singularities and to see the two elliptic curves emerge in d = 2, we reparametrise the problem through the Landau variable s = − (1−z) 2 z with s = (p 1 + p 2 ) 2 and we set m = 1 for simplicity. We then perform a rotation on the integrals I 4 , I 5 , I 6 , I 7 to disentangle the two elliptic curves such that their maximal cuts satisfy a GM system in block form The two blocks correspond to the elliptic Gauss-Manin differential system of the two-loop equal-mass sunrise evaluated at z and 1/z. For more details, in particular, how this system can be solved in d = 2, we refer to [66]. By simple ϵ-rescalings of the integrals I 1 , I 2 , I 3 and I 8 together with a normalisation by (1−z 2 )/z of I 3 , the corresponding differential equations are brought in ϵ-factorised form. To achieve an ϵ-factorised form also for the remaining master integrals I 4 , . . . , I 7 we apply our procedure. First, we need the semi-simple part of W corresponding to the Gauss-Manin system (5.4). Here we simply take the semi-simple parts of the two equal-mass sunrise Wronskians, is the holomorphic solution of the first sunrise block (c.f. eq.(3.27)) and is the other holomorphic solution of the second sunrise block. With this ingredient, we construct the rotation matrix to an ϵ-factorised basis. First, we take the inverse of W ss followed by a suitable ϵ-rescaling. Subsequently, all contributions of ϖ ′ 0 can be removed by another rotation. However, in contrast to the two-loop sunrise graph with equal masses, we need to introduce another rotation to factorise ϵ in the part of the differential equations where the two sunrise blocks couple to each other: where we introduced a new function defined by We expect that this function should be related to a new differential form defined on the geometry of the three-loop ice cone graph. Finally, the non-ϵ-factorised terms in the contributions from the subsectors can be removed by a final simple rotation. The final ϵ-factorised differential equation can be found in appendix B eq. (B.1).

The three-loop banana
It is by now well known that the geometry underlying the l-loop banana graph is an (l − 1)dimensional Calabi-Yau variety [40][41][42][43][44][45][46][47]88]. It is, therefore, interesting to see if our approach can also help to find ϵ-factorised bases when such more complicated geometries are involved. As the simplest case, we consider the equal-mass three-loop banana family (see figure  7) with propagators and numerators: A generic integral is then given by eq. (2.1). As customary when dealing with two-point functions, it is convenient to study the integrals in this family in d = 2 − 2ϵ dimensions and as functions of the single dimensionless parameter z = m 2 /p 2 . For the master integrals, we take the dotted basis given by I 1 = I 1,1,1,0 , I 2 = I 1,1,1,1 , I 3 = I 2,1,1,1 and I 4 = I 3,1,1,1 . (5.10) Notice that to simplify our notation we drop the indices corresponding to the additional numerators since we do not need them in the following. Moreover, we set m = 1 to shorten our formulas later. These master integrals satisfy a non-ϵ-factorised GM system d dz I = GM I , where the explicit form of GM is given in appendix B eq. (B.2). The maximal cuts of the three-loop banana master integrals satisfy a GM system which has as solutions the periods of a K3 surface [11,59]. As it is well known, differential operators of one-parameter K3 surfaces are symmetric squares of elliptic operators [89,90]. This was demonstrated explicitly for the differential operator of the three-loop banana graph in [11,88,91]. The homogeneous (3 × 3)-system admits the following solutions which can be, if desired, expressed through squares of elliptic integrals as explained in the aforementioned references. For our discussion, this particular form is not needed. Moreover, we hide the explicit dependence of z on the solutions to simplify our formulas. We can construct the Wronskian matrix Due to the underlying K3 geometry we get the unipotent part and τ -parameter From this we can construct the semi-simple part of W 3L given by (5. 16) We can use the Griffith's transversality conditions for the three-loop banana graph [80]: to eliminate in eq. (5.16) ϖ ′′ 0 in favour of ϖ 0 , ϖ 0 : Now we can continue as for the equal-mass two-loop sunrise family. We construct our rotation matrix starting with the inverse semi-simple piece of the Wronskian followed by an appropriate ϵ-rescaling. Notice, that by using the relation in eq. (5.18) we have already removed all contributions of ϖ ′′ 0 in the GM equations. We can then proceed to also remove all dependencies of ϖ ′ 0 in the GM system. So far, this is exactly the same procedure as in the two-loop case. Nevertheless, for the three-loop banana, this is not sufficient to achieve full ϵ-factorisation and we have to include an extra rotation where two new functions G 1 and G 2 have been introduced. They are defined as iterated integrals and bring the GM equations into ϵ-factorised form (see appendix B eq. (B.3)). This is exactly the same form as found in [50]. Also there the introduction of two new functions given by integrals over the holomorphic solution ϖ 0 was necessary. Notice that the new function G 1 only appears in the rotation matrix but not in the final GM matrix such that only expressions having simple poles appear in the ϵ-factorised GM matrix as also noticed in [50].

Conclusions and outlook
In this paper, we have elaborated on a procedure to obtain ϵ-factorised systems of differential equations for Feynman integrals which can be expressed in terms of iterated integrals defined on non-trivial geometries. We started by summarising physically and mathematically motivated criteria to select a good basis of master integrals from which to start our procedure. We have stressed that these criteria do not uniquely constrain the starting basis, but allow us to choose a set of integrals for which the subsequent steps of our procedure can be applied more easily. Once a starting basis has been identified, we proceed bottom-up, sector by sector. For each sector, the crucial step in our approach consists in constructing the Wronskian of the corresponding homogeneous differential equations and then rotating the basis by its semi-simple part. This construction was inspired by findings made by direct evaluation of elliptic Feynman integrals in terms of pure elliptic multiple polylogarithms [58]. After having rotated by the semi-simple part, one can then proceed by first cleaning up the homogeneous part of the system and then focusing on the inhomogeneous part. These clean-up steps can be seen as a generalisation of the standard approach used for polylogarithmic integrals and described in detail, for example, in [63]. The ϵ-factorised form that we achieve never contains derivatives of the periods, which we consider an important starting point to prove that the differential forms are independent, at least on the field of algebraic functions. An important feature of our approach is that, during the clean-up steps, one might have to introduce new non-trivial objects, which we refer to as G i and cannot simply be expressed as linear combinations of algebraic functions and of the periods of the geometry considered. These functions can instead be expressed themselves as iterated integrals over periods and algebraic functions and they can be related to similar objects identified in the differential equations for the multi-loop banana graph [51,52]. In our analysis, we have shown that, in many cases, the appearance of these new quantities can be related to either extra master integrals in the homogeneous system for ϵ ̸ = 0, or to extra singular points stemming from subtopologies.
We have successfully applied our procedure to various problems of increasing difficulty, starting from single-scale elliptic integrals, moving to multi-scale elliptic problems and also to families of integrals which involve geometries beyond an elliptic curve. In the elliptic case, the interpretation of our results is particularly transparent. As we verified explicitly in some of the examples above, we expect that a new function G i should be required for each independent differential form of the third kind that can be defined on the elliptic curve. This can either happen if there are extra master integrals in the top sector, which are linearly dependent in even integer numbers of dimensions, or if the extra residues are generated by the subtopologies. An interesting case in this respect is that of the non-planar three-point functions studied in section 4.2. Here, despite the presence of a third master integral in the top topology corresponding to an extra residue, no extra function G i is required. Indeed, we could prove that evaluating that differential form on one of the independent cycles, the integral reduces to a linear combination of the period and the quasi-period of the curve and it is therefore not linearly independent.
As already discussed, our procedure can also be applied to cases beyond elliptic. We considered as examples the three-loop ice-cone graph and the three-loop banana graph. Interestingly, also in those cases, the same steps take us to an ϵ-factorised basis and make it possible to identify in a very clean way a minimal number of new kernels that have to be added to the set of independent differential forms.
While many questions remain to be answered, we believe that our way of constructing ϵ-factorised bases of master integrals can not only provide a useful tool to solve practical problems but also help shed some light on the general properties of Feynman integrals beyond the polylogarithmic case and their evaluation in terms of iterated integrals over independent sets of differential forms. As the next steps, we hope to be able to study more in detail the extension of our procedure for either the Calabi-Yau case or higher genera geometries.
Here we provide the rotation matrix T sub which brings the polylogarithmic part of the second triangle family in canonical form The ϵ-factorised GM matrix of the second triangle family is given by -44 -

B Appendix C: Explicit results for ice cone and banana graphs
We give in this appendix some explicit results for the ice cone and banana graphs. The ϵ-factorised GM matrix for the three-loop ice cone is given by whereas the ϵ-factorised GM matrix is given by