On epsilon factorized differential equations for elliptic Feynman integrals

In this paper we develop and demonstrate a method to obtain epsilon factorized differential equations for elliptic Feynman integrals. This method works by choosing an integral basis with the property that the period matrix obtained by integrating the basis over a complete set of integration cycles is diagonal. The method is a generalization of a similar method known to work for polylogarithmic Feynman integrals. We demonstrate the method explicitly for a number of Feynman integral families with an elliptic highest sector.


Introduction
The method of differential equations [1][2][3] is the primary way phenomenologically relevant Feynman integrals are computed.This method works by taking the derivative with respect to a kinematic variable (or to several) of the members of a minimal basis of Feynman integrals for the problem, which we will refer to as master integrals.The result of taking that derivative may then be mapped back to the original master integrals using IBP identities [4,5], as implemented in for instance the computer codes FIRE [6], Kira [7], and others [8][9][10].This will result in a system of coupled linear differential equations for the master integrals, which may be solved with traditional methods.A breakthrough was made with the realization [11] that for a large class of Feynman integrals in dimensional regularization (with d = d 0 −2 ), the solution of the system gets much simpler, in some cases trivialized, by choosing a special basis with the property that the differential equation system is epsilon factorized.Denoting the basis of master integrals as J, this means that the differential equation system in the kinematic variable x may be written as where A (x) is a matrix dependent on the kinematics but independent of .
In the cases discussed in ref. [11], these matrices have the additional property that the individual A (x) -matrices may be unified as one matrix A (x) = ∂M/∂x with the property that the entries of M consist solely of logarithms of algebraic functions of the kinematics.A basis of master integrals for which the differential equations have these two properties (epsilon factorized, dlog form) is known as a canonical basis.
It is a well known fact that not all Feynman integrals can be brought to canonical form.One class (the simplest) of integrals for which a canonical form is unobtainable are those referred to as elliptic, a term that comes from the presence of elliptic curves at a cut surface.Such elliptic Feynman integrals, of which the fully massive sunrise integrals are a paradigmatic example, have been a subject of intense study in recent years .The purpose of this paper is to extend the epsilon factorized property of the differential equations given by eq.(1.1) to such elliptic cases.This has been achieved in the literature for a few specific examples [54,57] already, but this paper will present a more general algorithm for how to obtain such a basis.
For integrals for which a canonical form is obtainable, several approaches have been developed to finding it [61][62][63][64][65][66][67][68][69][70][71][72][73][74].One popular approach [61,67,68] includes the requirement that the master integrals are pure [75], which means that the value of the integral on each maximal cut (here in the sense of combinations of unitarity cuts that fix all degrees of freedom) is just a number, as opposed to a function of the kinematic variables.For elliptic Feynman integrals performing a maximal cut in that sense is impossible due to the presence of the elliptic curve.Yet a natural generalization of that cutting operation is an integral over a cycle, and such integrals can be performed also for the non-zero genus surfaces appearing for elliptic Feynman integrals (and beyond).
The algorithm proposed in this paper is the natural generalization of that maximal cut procedure, or more specifically it is to require that each master integral is non-vanishing on one and only one of the (members of a set of) basic cycles.This is similar to the prescriptive unitarity approach proposed in refs.[60,76] but here performed at the integral level for the fully dimensionally regulated integrals.
The use of d-dimensional unitarity cuts to investigate Feynman integrals and simplify their relations, has a long history.This was first done in the context of dimension shift relations [77], and later extended to differential equations [36], IBPs [78], and more [67,79,80].In some of those works [78,79] the notion of using a diagonal period matrix as a condition for obtaining an epsilon factorized differential equation, a notion central to this work, was prefigured, and so was the suggestion to keep certain epsilon-dependent prefactors before doing the cut analysis [67].Thus this paper may rightfully be regarded as the newest entry that series of works.
The dual vector space nature of the set of integrands ϕ i and integration contours γ j [78,[81][82][83], was clarified in recent work on the connection between Feynman integrals and the mathematical discipline of intersection theory [84][85][86][87][88][89][90][91], where the contours and integrands are viewed as representatives of a twisted de Rahm homology and cohomology respectively.This change of perspective makes it natural to look at the period matrix P ij ∼ γ j ϕ i of pairings between integrands and contours, and the algorithm proposed here is naturally formulated in terms of that object.
After introducing the notation in higher detail, we will in sec. 2 do two examples of the algorithm used in non-elliptic cases, while reformulating it in a way that generalizes to the elliptic case.Then in sec.3 we will make the complete formulation of our algorithm, and in sec.4 do a number of successful examples of its use for which many of the resulting expressions are added as an ancillary file.Finally in sec.5 we will discuss some loose ends and open questions, and summarize and conclude.In appendix A we will redo one of the examples with a different approach.

Notation and conventions
We will in this work encounter the complete elliptic integrals of respectively first, second and third kind.They are defined as where we use the Mathematica convention in which the arguments are k 2 and n 2 .The three complete elliptic integrals have the property that any integral of the form C R(x, Y )dx can be written in terms of them (along with elementary functions), where R is any rational function, Y is the square root of a polynomial in x of degree three or four, and where C is any closed contour.We will in this work only discuss two-loop Feynman integrals, which can be written as where the propagators D i are quadratic functions of the momenta on the form k 2 −m 2 and d is the spacetime dimension.As mentioned earlier, we will consider the integrals in dimensional regularization where d = d 0 − 2 with d 0 being an integer, in this work either 2 or 4.
We will consider our Feynman integrals in Baikov representation, in either its standard [92][93][94] or loop-by-loop [67,95] version.Thus the integral may be written where the x i are Baikov variables (that equal the original propagators, x i = D i ), the α j are generally irrational powers, and where u is a multivalued function that will have the general form where the B j are polynomial functions called Baikov polynomials of the x i .The prefactor K will generally be given in terms of gamma functions of arguments linear in the spacetime dimension.In the examples we will look at, K will be a pure function, by which we will mean a function with the property that the coefficients in the expansion has uniform weight [96] such that the coefficient of n has weight one higher than that of the coefficient of n−1 and only numerical prefactors.We will not make any attempts at generalizing these notions to elliptic objects.
We will in this paper make significant use of generalized unitarity cuts, in particular maximal cuts which refers to cuts of all genuine propagators (and thus differs from the notation in the N =4 literature, in which a maximal is defined to fix all degrees of freedom).This is due to the fact that differential equations [36], as well as IBP relations [78], remain the same after such a cutting procedure, and indeed the Baikov parametrization is particularly suited for generalized unitarity cuts at the integral level [67,80].
In the examples in this paper, the maximal cut will leave a single integration to be done, bringing the integrals to a univariate form of eq.(1.6) where z is the remaining Baikov variable, and φ is a rational function of z and the kinematical variables.We will occasionally use the abbreviation φ := φdz.In general this notation follows that of the recent work on intersection theory and Feynman integrals such as refs.[84,85].We will also discuss various integration contours.We will name them with the symbols C or γ, and use the notation C i for a contour surrounding a pole at z = i, and the notation C i-j for a contour surrounding a branch cut between the points z = i and z = j.

Notation -The sunrise integrals
We will in this paper work quite a lot with the sunrise integral with different mass distributions.Rather than repeating the same definitions over and over, let us write them here in the most general case: The three-mass sunrise integral is defined by the complete set of propagators with p 2 = s, and where only the first three may appear as genuine propagators.With these definition, the integrals are given as We will always look at the case where d = 2 − 2 , when discussing the sunrise.
In the loop-by-loop Baikov representation, this integral may be written as where is the product of Baikov polynomials raised to their corresponding powers, and where is a pure function.The choice a 5 = 0 made in eq.(1.11) is done in order to simplify the loop-by-loop parametrization, which in this case is performed by first integrating out the k 2 -loop.This limitation of the loop-by-loop Baikov parametrization can be overcome (for instance using a Passarino-Veltman [98] inspired approach as described in section 11.1 of ref. [85], and for a recent alternative approach to this see ref. [99],) but that additional complication is not necessary for the present discussion.On the maximal cut (i.e. of D 1 , D 2 , D 3 , and with z = D 4 ), u reduces to meaning that any integral in the sunrise family can be written as on that cut, where φ will be some rational function of z, the exact form of which will depend on the values of the propagator powers a i .We note that u| 3×cut is a pure function divided by the square root of a polynomial of degree four, a fact we will be using implicitly in the following.

Motivation
In this section we will discuss two simple examples of non-elliptic Feynman integrals, and show how to find their canonical forms with the traditional method.This method will then be reformulated in terms of diagonalization of the period matrix, giving a formulation that generalizes directly to the elliptic case.

The double box
Let us as an illustration of the procedure discuss the well-known fully massless double box.That family of Feynman integrals (first computed in [100,101]) was used as the example in ref. [11] which introduced the concept of epsilon-factorized differential equations.It is given by the set of propagators of which only the first seven are allowed to appear as genuine propagators.
The family contains eight master integrals of which two are in the highest sector, and it is those two that will interest us in the following.In ref. [11] it was found that a set of integrals giving epsilon-factorized differential equations are where c = 4 s −2 is a pure prefactor to be ignored in the following.How can we know a priori the two prefactors s 2 t and s 2 ?One way is to perform the maximal cut (in the sense of cutting the seven actual propagators) of the integrals in the family.Using the loop-by-loop Baikov parametrization, this gives where is a pure function.
The prescription is now to factorize out the pure part of the integrand, which here corresponds to taking the limit → 0 of the object under the integral sign.This (ignoring the already pure K) gives the integrand Continuing the cutting procedure shows us that I db 111111100 has two residues.One at z = 0 evaluating to 1/(s 2 t) and the other at z = t evaluating to −1/(s 2 t).On the other hand

I db
1111111 −10 also has two residues, one at z = t evaluating to −1/(s 2 ) and the other at z = ∞ evaluating to 1/(s 2 ).So we see that the two integrals become pure if I db 111111100 is given a prefactor of s 2 t and I db 1111111 −10 a prefactor of s 2 , as it is given by eq.(2.2).A different way of formulating the procedure is by writing up the period matrix for the two integrals, again after factorizing out and discarding the pure part of u.We pick the two integration cycles as γ 1 being the cycle C 0 surrounding the (z=0)-pole and γ 2 being C ∞ surrounding the (z=∞)-pole, with the various contours shown on fig. 2. Giving the two integrals prefactors f 1 and f 2 we get the period matrix and requiring the period matrix to be 2πi times the unit matrix I fixes the f i to the values found above.
Making this procedure more general and algorithmic, we could also write each of the canonical integrals we want to find as general linear combinations of the two precanonical intermediate basis integrals I db 111111100 and I db 1111111 −10 , that is This gives the period matrix and once again imposing P = 2πiI fixes the f ij coefficients uniquely to again corresponding to the prefactors given by eq.(2.2).

The two-mass non-elliptic sunrise
Let us look at another non-elliptic example that looks more similar to the elliptic examples we will encounter later.That is the two-mass non-elliptic sunrise integral (sne), defined as in eqs.(1.9) and (1.10) but with m 2 3 = 0 and m 2 1 = m 2 2 ≡ m 2 .Using the loop-by-loop Baikov parametrization of this integral as in eq.(1.11), we may make a univariate representation of the integral on the maximal cut.Following our If one were to follow the prescription of cutting the last variable in order to impose a maximum cut of ±1, one would (formally) first have to get rid of the square root through rationalization.This can be done by a variable change such as but let us try to proceed without such procedures.The integration plane can be seen depicted on fig. 3.One can define a basis of independent cycles as a cycle γ 1 surrounding the (z = s)-pole, and one γ 2 surrounding the pole at infinity.Picking the contour around the branch-cut from z = 0 to z = 4m 2 as one of the master contours would be equally valid, but the above choice gives the nicest result.We will make a choice of intermediate basis integrals (precanonicals) as I sne 11100 and I sne 111 −10 .This corresponds to the integrands φ1 = 1 and φ2 = z (2.12) With this we may compute the integrals on the cycles, giving Writing the two integrals we are looking for as generic linear combinations we obtain the period matrix P ij = f ik g kj as Solving for P = 2πiI gives corresponding to (2.17) With this we may compute the epsilon factorized differential equation matrix written disregarding integrals in lower sectors (here only the double-tadpole I sne 11000 ) so we see that our method works in this case as well.This success is what motivates us to try the same approach in elliptic cases.

The proposed algorithm
Motivated by the examples in the previous section, we are now ready to formulate our algorithm.Having a set of ν integrals of the form where as discussed in section 1.1 u is a multivalued function, and φi a set of ν rational (in x) functions, our claim is that if the integrand can be written as where σ is a pure function, and where the period matrix P with P ij = γ j Φi d n x is proportional to the identity matrix as P = (2πi) n I, then the integrals J i will fulfill epsilonfactorized differential equations.In cases where K is itself pure to begin with, it can be ignored in the above discussion.What we will do in practice is to write the basis integrals as a linear combination of known intermediate basis integrals I j (often called precanonical in the polylogarithmic case,) that is J i = f ij I j or correspondingly as Φi = f ij Φint j , and then imposing P = (2πi) n I will impose ν 2 constraints on the f ij fixing them all uniquely.
In the examples in the following section, we will write the Φint as Φint i = φi /Y where Y is the square root of a (monic) polynomial of degree four, corresponding to the elliptic curve characterizing the problem.

Examples
In this section we will look at examples of elliptic Feynman integrals which we can bring into a form that has epsilon-factorized differential equations with the above algorithm.We will restrict the discussion to cases where the ellipticity is present in the highest sector only, and where a maximal cut can bring the integrals in that sector to a univariate form.

The same-mass elliptic sunrise
As a first example let us look at the most basic of elliptic Feynman integrals, the same-mass elliptic sunrise (s1m).This is defined as in section 1.2, but with the restriction that the three masses are the same, i.e. m 2 1 = m 2 2 = m 2 3 = m 2 .This integral family has three master integrals.One is the double-tadpole (e.g.I s1m 11000 ) while the last two are in the highest sector and the only ones we will care about in the following where we disregard subsectors completely.Picking as intermediate basis integrals 11100 and I s1m 21100 we can write the two candidate integrals as where the "lower" covers potential contributions from the double-tadpole sector.The c is a constant prefactor, the value of which will not change the differential equations and which will be ignored in all of the following.
Our goal now is to find values of the f ij that makes the period matrix of these two integrals the identity matrix.
Parametrizing this integral with the loop-by-loop Baikov parametrization, we get from eq. (1.11) as by eq.(1.12).K given by eq.(1.13) is a pure function, and u is pure multiplied with a factor of Y −1 where on the maximal cut (i.e. of x 1 ,x 2 ,x 3 and with Following the algorithm of sec. 3 we factor out the pure part of u, leaving the expressions where the integrands are given by of which the latter has been computed as φ2 = (∂ x 1 u)/u| cut .Please note the presence of the -dependence in φ2 .This factor comes out of the algorithm in a natural way, and is essential in order to get the correct epsilon factorization Y 2 has four roots so we can look at the two basic cycles γ 1 = C ii-iii and γ 2 = C i-ii which may be seen on fig. 4. In the convergent case allowing us to compute the integrals and The period matrix (  Imposing M = 2πiI gives a unique solution for the f ij as With this choice we can compute the system of differential equations for the two integrals.It is epsilon factorized with and a similar epsilon-factorized differential equation may be found in the other variable dJ/dm 2 .This shows that our algorithm works for elliptic examples too.

The elliptic nonplanar double triangle
The elliptic nonplanar double triangle (npt) [41,102], is defined by the set of propagators where only the first six may appear as actual propagators.The kinematics is such that p2 1 = p 2 2 = 0 and (p 1 +p 2 ) 2 = s.This sector contains 11 master integrals.The first nine are in subsectors, and can be put in canonical form in the traditional sense.The last two are in the highest, elliptic sector, and it is those that will concern us here.Choosing as intermediate basis integrals the set I npt 1111110 and I npt 2111110 we will write the two elliptic master integrals as linear combination of the intermediate basis: Using the Baikov parametrization, the integrals in this sector may on the maximal cut (of x 1 -x 6 with z = x 7 ) be written as Around d = 4, we have once again that K is pure, and u is pure times a factor of Y −1 with where a factor of s has been absorbed in the φ.This justifies studying the integrals with the u → Y −1 limit taken, and in that case we may write the intermediate basis integrals as where The four roots of Y 2 are given as and we see that a set of independent integration contours can be chosen as similarly to the previous example where the contours are depicted on fig. 4. We note that where in the previous section we always had C α-β φ = 2 β α φ, that is not the case here for C i-ii φ 2 /Y due to a divergence in z = r i .The results for the integrals are where we have We may now compute the period matrix P ij = f ik g kj , and imposing P = 2πiI allows us to uniquely fix the coefficients f ik of eq.(4.17) to With this we may compute the differential equation, which once again is epsilon factorized and it is given by with with a similar and likewise epsilon factorized equation to be found for the other derivative dJ/dm 2 .

The two-mass elliptic sunrise
Next let us look at the two-mass elliptic sunrise (s2m).This not very well studied integral is defined as is section 1.2 with m 2 1 = m 2 2 = m 2 b and m 2 3 = m 2 a , making it intermediate between the same-mass elliptic sunrise of section 4.1 and the three-mass elliptic sunrise of section 4.4 both in terms of the number of scales and the number of master integrals.The integral family contains five master integrals; two double-tadpoles and three in the highest elliptic sector which is what will concern us here.
On the maximal cut in d = 2 we get that the integrals in the highest sector may be written as We will pick the three intermediate basis integrals I s2m 11100 , I s2m 21100 , and I s2m 111 −10 , which correspond to the integrands Once again we will define the master integrals we are looking for, as generic linear combinations of the intermediate basis integrals Y 2 has the four roots and with this we may define our three master contours as depicted on fig.6.We will also define the abbreviations and where we note the relation With this we can compute the needed integrals, as g ij = γ j φi dz Y : making this the first appearance of the complete elliptic integral of the third kind Π(n 2 , k 2 ).With this we may write down the period matrix P ij = f ik g kj and impose it to be diagonal P = 2πiI by solving for the f ik .We get δ + (1+2 ) (4.37) Once again the differential equation matrices are epsilon factorized dJ i /dx = A (x) ij J j .The entries are in general too large to be written here, one example entry is and for the full expressions see the added file.

The three-mass elliptic sunrise
Next we will look at the generic or three-mass elliptic sunrise (s3m), exactly as it is defined in section 1.2.The integral family contains seven master integrals; three double-tadpoles and four in the highest, elliptic sector.Integrals in that highest sector may, on the maximal cut and in two dimensions, be written as We will pick our master integrals as linear combinations of now four intermediate basis integrals, as These four intermediate basis integrals correspond to the integrands where the latter, which may appear incompatible with the loop-by-loop parametrization used, has been obtained using the Passarino-Veltman inspired approach described in section 11.1 of ref. [85]. and where again we have the relation (1−k 2 )/n 2 = 1 − k 2 /n 2 and likewise for ñ.With this we may write the integrals h ij = γ j ξ i /Y : We may then form the integrals of the φ integrands g ij = γ j φ i /Y as combinations of the above, and then form the period matrix P ij = f ik g kj .Fixing it to be diagonal P = 2πiI will fix all 16 f ij .We will not list them all here, two examples are With this choice we get the differential equations in epsilon-factorized form dJ i /dx = A (x) ij J j for all the kinematic variables x.The matrix entries are too large to be written here.

The added expressions
In many of the examples done in this section, the result were deemed too large to be written out in full.For that reason a file expressions.mhas been added to the arXiv version of this paper, which contains all the expressions in Mathematica format.The expressions are named for instance Ass2m for the matrix A (s) s2m or frulesnpt for replacement rules for the f ij in the case of the non-planar double triangle.To obtain the full list of expressions defined in the file one might use the Mathematica command Names["Global'*"] after reading it in.

Form and integration of the differential equations
One topic that has been avoided so far in this paper, is how to solve the differential equations once the epsilon-factorized form has been obtained.For Feynman integrals in the traditional canonical form it is often straight forward to integrate them up order by order in the epsilon expansion, giving results in the function class of generalized polylogarithms [103], (at least in the case where all square-roots can be rationalized, which is however known to not always be possible [104]).Many attempts have been made to extend this function class to elliptic cases and beyond [22,28,40,46], yet none of those seem directly applicable to the form of the differential equations derived in this paper.Further investigations will be needed into how to best proceed analytically from the differential equations found here.
If, however, one is satisfied with a numerical approach, the development of techniques to numerically integrate this type of differential equations is developing quickly [49,[105][106][107][108][109][110].Not all of these approaches require the differential equations to be epsilon factorized, but they all benefit from from it and work faster in that case, so even in the worst case scenario the techniques developed here can be used to speed up the numerics.
We do however notice one property that the differential equations found on the previous pages share with those in traditional canonical form and which might cause optimism for the prospect of analytical integration: They are free of higher poles at singular points including the point at infinity.This property is not obvious, it might appear as if for instance the A 11element of the differential equation matrix found for the same-mass sunrise of eqs.( 4 and we see that no higher pole in y is present, the candidate having canceled with a term from the expansion of the elliptic integrals.

Freedom in basis choice
There is some amount of freedom in the algorithm discussed on the previous pages.This freedom corresponds to performing row operation on the period matrix.For instance one might define the master contours γ i as some linear combination of the ones used in this paper.In all the elliptic examples we chose γ 1 = C ii-iii and γ 2 = C i-ii , but one might equally well have chosen for instance γ 2 = C iii-iv which would yield a slightly different set of integrals in the end, and also there is nothing preventing the choice of some C-linear combination of contours, such as γ = 5 7 C ii-iii + 8i 13 C i-ii .Yet these kinds of redefinitions seem to provide no simplification of the integrals or their differential equations.Likewise one might take linear combinations of the resulting integrals, since if the set J has epsilon factorized differential equations so will a C-linear combination J = BJ with B ij ∈ C. For example for the case of the two-mass elliptic sunrise, it seems from eqs. (4.37) as if it would provide a slight simplification to perform the replacement J s2m 2 .Yet as mentioned in the introduction, there are previous examples in the literature in which elliptic Feynman integrals have been put into a form that allows for epsilon factorized differential equations, and which are not equal, or in a form related through the above mentioned transformations, to the expressions found here.In particular the same-mass sunrise of section 4.1 has been discussed in ref. [54], and the three-mass sunrise of section 4.4 has been discussed in ref. [57].It must be admitted that the integrals found in those references are "nicer" than those found with the algorithm outlined here.In particular the matrix of f ij coefficients are in both cases found to be something more akin to a lower triangular form and the differential equation matrix is nicer as well as we will see.Focusing on the same-mass sunrise, the form found in ref. [54] has f 12 = 0 which means that the master integrals there may be written where the values are found to with the definitions of k and R being as given in sec.4.1.We see that these expressions are quite different from what was found in eqs.(4.13).Aside from having three coefficients rather than four, eqs.(5.3) has simpler kinematics dependence, and only complete elliptic integrals of k 2 appear where eqs.(4.13) also contained complete elliptic integrals of argument (1−k 2 ).On the other hand the π and dependence is simpler in eqs.(4.13).
From this we may compute the differential equations.They may be written as d J/ds = A Js1m with which is simpler than the expressions in eqs.(4.15) by a substantial amount, we see for instance that the only elliptic integral appearing is K(k 2 ), and that the diagonal entries are free of even this.On the other hand the K(k 2 ) appears both in numerator and denominator, where in eqs.(4.15) the elliptic integrals appeared in the numerator only.
We may compute the period matrix of this example.The result is P ij with and in addition it is worth noticing that the determinant of this matrix is given as We see that P 11 and det(P ) both are constants, in the sense that they are free of any kinematic dependence including through elliptic integrals.This may be used as a guideline for obtaining an -factorized form alternative to the one discussed in this paper.We do an example of this in Appendix A, but it is not clear if this can be generalized to elliptic sectors with more than two master integrals.
It would be great if some principle could be found to a priori generate the transformation between the forms discussed in this paper, and forms similar to the one discussed above.

The number of cycles
In principle there is a one-to-one correspondence between the number of master integrals and the number of independent cycles [82,83,86].This correspondence is utilized by the Lee-Pomeransky criterion [82], which states that the number of master integrals ν is given as ν = number of solutions to "ω = 0" where ω := dlog(u) (5.7) We can test this for the examples discussed in secs. 2 and 4 of this paper.The result is that the numbers agree in all the cases, except for the same-mass sunrise of sec.4.1 and the nonplanar double triangle of sec.4.2, which are the two cases in which an elliptic sector has two master integrals.For the case of the same-mass sunrise, this "miscounting" is well known, and is discussed in further detail in refs.[85,89] where it is shown to be caused by peculiarities in the interplay between the loop-by-loop Baikov parametrization and the maximal cut.Focusing on the same-mass sunrise it might seem as if a valid third master integral would be I s1m 111 −10 corresponding to φ3 =z, with a corresponding third master contour surrounding the pole at infinity γ 3 =C ∞ making the set of cycles look as in the case of two-mass elliptic sunrise of fig. 6.Yet we know from IBP relations that and it might be instructive to see how this relation is realized over the basic cycles.Using the φi and γ j from sec.4.1, and again defining g ij = γ j φi dz/Y , we get ).On the other hand for contour γ 3 the relation clearly does not hold since g 13 = 0, while g 33 = −2πi which is then the amount with which the relation is broken.At the last contour γ 2 we get using a similar relation for a special point of Π(n 2 , k 2 ), that g 32 − 3m 2 +s 3 g 12 = 2πi 3 (5.12)so we see that when the peculiarities of the loop-by-loop Baikov parametrization ruins the counting of master contours, relations such as eq.(5.8) break only by a factor proportional to 2πi.
For the nonplanar double triangle of sec.4.2 the behaviour is extremely similar.There the extra IBP-derived integral relation relating an integral with a pole at infinity to the original two, is −s 2 I npt 1111110 + lower (5.13) and that relation holds exactly on contour γ 1 and is broken with factors proportional to 2πi on γ 2 and a potential γ 3 surrounding a pole at infinity.The main difference is that the Π(n 2 , k 2 ) special value relation that makes this explicit is even simpler in that case.

Freedom in intermediate basis choice
One might consider if the intermediate bases chosen in the examples discussed in secs. 2 and 4 have some special property that allows for the algorithm to go through that would not be shared with any arbitrary basis choice.The answer to this is mostly no, any set of valid master integrals might be chosen as the intermediate basis.Let us illustrate that with an example, that of the two-mass elliptic sunrise of sec.There are however cases in which not any choice of intermediate basis is equally convenient.This is true in the cases that have the "miscounting" of independent cycles discussed in sec.5.3, such as the same-mass elliptic sunrise.There we saw that a counting of independent cycles seems to include a third master contour surrounding a pole at infinity such generalizations in this work, but I hope that my approach contributes with insights that may help bring us to a point where such a step might be taken.
The algorithm presented in this paper is able to systematically bring differential equations for elliptic Feynman integrals into a form for which the differential equations are epsilon factorized.We demonstrated this for a number of examples of varying complexity, and it is the hope that this development will be one step on the way towards bringing the understanding of elliptic Feynman integrals to the same level as the polylogarithmic case.

Figure 2 .
Figure 2. Contours discussed for the double box

Figure 3 .
Figure 3. Contours discussed for the two-mass non-elliptic sunrise

Figure 4 .
Figure 4.The two independent contours for the same-mass elliptic sunrise

Figure 6 .
Figure 6.The three independent contours for the two-mass elliptic sunrise

Figure 7 .
Figure 7.The four independent contours for the three-mass elliptic sunrise

4 . 3 . 20 ( 5 . 14 )
Instead of the choice of (I s2m 11100 , I s2m 21100 , I s2m 111 −10 ) made in sec.4.3, one might instead pick for instance the set (I s2m 11100 , I s2m 111 −10 , I s2m 111 −20 ) which is the choice naturally made by Kira.With this one can go through the same procedure as before, define Ji = fi1 I s2m 11100 + fi2 I s2m 111 −10 + fi3 I s2m 111 −compute the period matrix, and impose it to be proportional to I thereby fixing the fij .The resulting Ji equal those from sec.4.3 up to integrals in lower sectors as well as terms of order :Ji = J i + O( ) + lower (5.15)Additionally the new integrals have epsilon factorized differential equations d Ji /dx = Ã(x) ij Jj in each of the three kinematic variables x ∈ {m 2 a , m 2 b , s}.The expressions for f and the Ã(x) can be found in the file added to the arXiv version of this paper.