An algorithmic approach to finding canonical differential equations for elliptic Feynman integrals

In recent years, differential equations have become the method of choice to compute multi-loop Feynman integrals. Whenever they can be cast into canonical form, their solution in terms of special functions is straightforward. Recently, progress has been made in understanding the precise canonical form for Feynman integrals involving elliptic polylogarithms. In this article, we make use of an algorithmic approach that proves powerful to find canonical forms for these cases. To illustrate the method, we reproduce several known canonical forms from the literature and present examples where a canonical form is deduced for the first time. Together with this article, we also release an update for INITIAL, a publicly available Mathematica implementation of the algorithm.


Introduction
The computation of Feynman integrals is of immense importance in perturbative Quantum Field Theory.Not surprisingly, a number of different approaches have been developed and polished to a high degree to tackle this difficult challenge.Among these methods, differential equations have proven to be one of the most powerful and therefore most widely used approaches [1][2][3][4][5].In this method, one first uses integration-by-parts identities to reduce a given set of Feynman integrals to a basis of so-called master integrals ⃗ f .Then, one computes the derivative of the master integrals with respect to the kinematic invariants.The result is a linear combination of Feynman integrals which, due to the integration-byparts identities, can again be written in terms of the basis ⃗ f .For one kinematic variable x, the differential equations therefore take the form where the coefficient matrix A(x, ϵ) is a rational function in x and the parameter of dimensional regularization ϵ.
Although eq.(1.1), together with a boundary condition, fully determines ⃗ f , the matrix A(x, ϵ) is often complicated, and it is difficult to solve the differential equations analytically.(For numerical approaches, see e.g.[6]).A paradigm shift occurred with the realisation that the conjecture [7,8] that (1.1) can be transformed into much simpler form, with the help of a basis transformation, ⃗ g = T ⃗ f , for a suitable invertible matrix T (x, ϵ).The two main features of (1.2) are that firstly, its RHS is proportional to ϵ, and secondly, the singularity structure of Ã(x) manifests the Fuchsian property of the system.The first feature means that the solution to (1.2) in terms of a power series in ϵ is reduced to straightforward iterated integration.The second feature restricts the class of iterated integrals, and in practice often allows one to fix the boundary constants in a simple way, see e.g.[9].For example, in the case of Feynman integrals evaluating to multiple polylogarithms (MPLs) [10,11]-a very important class of special functions in this area of research-the canonical differential equations (1.2) take the specific form for some set of rational functions α i (x) and constant (kinematic-and ϵ-independent) matrices m i .Indeed, the form (1.2) together with (1.3) makes it manifest that resulting iterated integrals are multiple polylogarithms (MPLs).
It is important to mention that the canonical differential equations can conjecturally be obtained from an analysis of the loop integrand of Feynman integrals [7,[12][13][14], prior to integration.This is closely connected to the property of uniform transcendental weight (see [15] and references therein).If the master integrals are chosen according to this integrand analysis, one immediately finds differential equations in a canonical form, without the need for constructing a (possibly complicated) transformation matrix T .See refs.[9,[16][17][18] for state-of-the-art applications.Alternatively, one may first compute A, and then try to algorithmically construct T .See refs.[8,9,[19][20][21][22][23][24] for various ideas in that direction, including powerful algorithmic implementations.The idea of a canonical form and concrete ways of obtaining it has streamlined the computation of Feynman integrals and led to significant advances in the computation of Feynman integrals, and corresponding physical applications in the last ten years.
However, with increasing loop order or increasing number of kinematic variables, there are many cases of Feynman integrals evaluate to functions beyond MPLs, and therefore eq.(1.3) needs to be generalized.The simplest of those cases is when the matrix Ã(x) involves a single type of elliptic integral satisfying a second-order differential equation.The natural question to ask is what the precise form of the canonical differential equations is in the elliptic case, and beyond.This question has received significant attention in recent years.
In the literature, several cases of differential equations in the form (1.2) can be found [25][26][27][28].These results are sometimes referred to as ϵ-factorized forms.This terminology is certainly correct in view of eq.(1.2).In this paper, we use both the expression canonical form and ϵ-form, almost interchangeably.For us, the term canonical form means that we look for a simplified system of differential equations as first conjectured in [7], while ϵ-form leaves open the possibility that further simplification to A(x) may be found in the future.Indeed, settling conclusively the question of the specific form of A(x) is a key open problem in this research area.
This important question is intimately linked to the class of iterated integrals one expects to appear in the answer.For example, in the case of MPLs, it is clear that the integration kernels in eq.(1.3) are sufficiently general to cover all cases.What generalizes these integration kernels for elliptic cases?In the literature for elliptic multiple polylogarithms, one finds a host of different representations.A popular class of iterated integrals involves modular forms (or various generalizations thereof) as integration kernels, see e.g.[29][30][31][32].The different orders in ϵ of the corresponding basis ⃗ g are then written in terms of iterated integrals of modular forms.For certain applications and integration contours, the latter class of functions is equivalent to the elliptic multiple polylogarithms (eMPLs) [30,[33][34][35][36], which are often used for the direct integration of elliptic Feynman integrals from their parametric representation.(Because of the importance of these two types of iterated integrals, the study of their properties, relations, analytic continuation and numerical evaluation is rapidly progressing [32,37,38].)However, despite these advances, a final picture for what integration kernels are needed in the canonical differential equations has not yet been established.
One approach to get insights into this is to extend the integrand analysis to the elliptic case.This has been pursued in [30,[39][40][41].In elliptic cases, after taking a certain number of residues, one encounters an elliptic curve.This implies that 'leading singularities' are now not just maximal residues, but rather correspond to independent integration cycles on that elliptic curve.For example, in the case of the sunrise integral, there are two independent integration cycles, which is why the key part of the differential equation is given by a coupled two-by-two system.The maximal cut integral by definition solve the corresponding secondorder equation, but going back to a first-order system leaves some freedom, i.e. additional information is needed to fix the canonical form.
In this paper, we follow a complementary approach.Leveraging information from maximal cuts of elliptic integrals, and imposing desirable properties such as Fuchsian behavior, we make an ansatz for what integration kernels the differential equations may contain, and then determine algorithmically whether a transformation to such a canonical form exists.Our goal is to extend the algorithm of [23] from the polylogarithmic case to the elliptic case.To determine the unknown coefficients, our algorithm assumes that at least one of the integrals is already pure, i.e. it should already be part of the canonical basis ⃗ g.We then require that the Picard-Fuchs equation derived from ⃗ f should equal the one derived from ⃗ g, which gives enough constraints to determine the ϵ-form and therefore the rest of the integrals in ⃗ g.
The paper is organized as follows: Section 2 is split into two parts.In section 2.1 we review the algorithm of [23] to set notation and clarify important concepts.Section 2.2 discusses how the ansatz has to be adapted to the case of elliptic Feynman integrals.Specifically, our approach will be that the integration kernels should have similar properties to the ones defined in [25,29].For this, we also take inspiration from already known ϵ-forms [25][26][27][28].We then provide several non-trivial examples in section 3.In the first example, we reproduce the known ϵ-form of the kite integral family [25].The goal of this section is to clarify how to extract the elliptic functions from a given set of differential equations.Then, in sections 3.2 and 3.3, we present previously unknown ϵ-forms for an example involving square-roots and linearly-dependent derivatives, respectively.In sections 3.4 and 3.5 we encounter examples where we have to adapt our ansatz to functions satisfying a third-order differential equation.Lastly, we describe our implementation in section 4 and conclude in section 5. Three appendices provide additional details to the material covered in the main text.

Description of the method
In this section, we first review the algorithm of [23], which itself is based on [42], and then describe how it naturally extends to the elliptic case.

Review of the algorithm
Given a basis of n master integrals ⃗ f = (f 1 , . . ., f n ) T depending only on a single scale x, one can use integration-by-parts identities to derive the system of differential equations The key idea of the algorithm is to assume that there is one integral, e.g.f 1 , which is already part of the canonical basis ⃗ g and does not require any further transformations, i.e. f 1 = g 1 .The first step is then to remove all dependence on the remaining integrals by transforming to the basis formed only by f 1 and its derivatives1 : 2) The transformation matrix Ψ is given by where ) and ⃗ v 1 = (1, 0, . . ., 0).In a second step, one can invert Ψ and project eq. ( 2.2) on the first row to obtain an n-th order differential equation for f 1 (Picard-Fuchs equation): where the coefficients are Next, one repeats the above steps with the matrix A(x, ϵ) replaced by B(x, ϵ), where the latter represents our ansatz for the canonical form: and the m i are constant matrices that will be determined in the following.In analogy with eq.(2.2), this process involves the basis change where the matrix Φ(x, ϵ) can now similarly be used to derive an n-th order differential equation for g 1 with coefficients given by −⃗ v 1 Φ −1 , c.f. eq.(2.7).From the assumption which can be solved for the unknown constant matrices m i .The transformation ⃗ f = T⃗ g is then given by T = Ψ −1 Φ.

Ansatz for the canonical form
Besides choosing a suitable initial integral, a key ingredient of our algorithm is the ansatz for the canonical differential equations (2.8).In the particular case where the result can be written in terms of multiple polylogarithms (MPLs), the integration kernels take the form where the α i (x) are rational (algebraic) functions that can often be determined through the poles of A(x, ϵ) (at least in the univariate case).Note that (2.11) has only simple poles.
It is expected that one can always find a form of Ã(x) such with at most a single pole at a given singular point.This corresponds to the Fuchsian property of Feynman integrals.For more information, see the discussion in refs.[8,19].In principle, one could relax this condition and allow for double poles and even transcendental functions to appear in the integration kernels, but this is not desirable.In particular, ideally the integration kernels directly lead to the desired class of iterated integrals, without further manipulation (like integration by parts or shuffle algebra).
Starting at two-loop level, there are many examples of Feynman integrals that evaluate to special functions beyond MPLs.To understand when this is the case, following [8], let us imagine we have found a basis ⃗ f in which the matrix A(x, ϵ) is analytic in ϵ, i.e.
and then remove the ϵ 0 -term The class of functions appearing in the solution is then expected to be a superset of the class of functions appearing in ϵ-form.For example, in the polylogarithmic case, the transformation T (0) (x) involves only rational (or possibly algebraic) functions in x, as well as MPLs, from which we then conclude that (2.11) is expected to be sufficient for the ansatz.We note that it is usually enough to apply the procedure of integrating out the ϵ 0 -part only to the diagonal blocks of the differential equations such that A (0) (x) becomes strictly lower triangular after the transformation T (0) (x).Our terminology will therefore be that a certain diagonal block (also called sector in the following) introduces a certain set of functions into the ansatz.After MPLs, the next most complicated case arises when the T (0) (x) for a specific sector involves a function satisfying a second-order differential equation where α 0 (x) and α 1 (x) are rational functions, see section 3 for explicit examples.In principle, our ansatz then has to include all linearly independent functions which are rational in x, Ψ 1 , Ψ 2 and the derivatives Ψ ′ 1 and Ψ ′ 2 : where the powers k j are, a priory, arbitrary integers and c l are the same singular points as in the differential equations matrix A(x, ϵ).To restrict this further, we again require that the integration kernels directly lead to a well-understood class of iterated integrals, namely the ones discussed in [29], see also [31].Note that this assumption has been observed to be correct for a large class of Feynman integrals involving complete elliptic integrals, however, in general we do not expect it to hold for all such Feynman integrals.In particular, we will discuss an example violating this assumption in section 3.5.Imposing this restriction leads to the following conditions on the integration kernels: • No double poles are allowed.
• Only one of the two solutions to (2.14) can appear, not both at the same time.For concreteness, we call this solution Ψ 1 .
• Only Ψ 1 and not its derivative can appear.2 • The minimum degree of Ψ 1 is −2.
While the first restriction reflects the well-understood Fuchsian property of the system, the remaining three conditions are imposed so that the integration kernels transform in a specific way under modular transformations, see e.g.[29,31].
Our ansatz is therefore subject to the condition that there are no double poles and we only take linearly independent combinations.Further, we expect that the maximum degree of Ψ 1 depends on the specific example and therefore we will come back to this point in the next section.We note that, for the purpose of our algorithm, it is not required to have an explicit expression for the function Ψ 1 , but it suffices to know the second-order differential equation it satisfies and to study its behavior near all singular limits such that one can restrict the ansatz to be free of double poles.This behavior can conveniently be obtained by either applying the method of Frobenius on (2.14) or the method of Wasow [44] on (2.13) (see also [45]).Lastly, we mention that the restrictions discussed in this section are sufficient for finding an ϵ-form of nearly all examples discussed in this paper after adjusting the singular points and the function Ψ 1 accordingly (with the banana integral family of section 3.5 being the only exception).

Examples and applications
As discussed in the introduction, we show the application of our algorithm to several nontrivial examples, where each of them introduces a specific new concept.As a warm-up, we first reproduce the known ϵ-form of the kite integral family.Then we discuss two new examples involving square-roots and linearly dependent derivatives, respectively.Lastly, we show on two examples the applicability of our algorithm to functions satisfying a third-order Picard-Fuchs equation.The first application we want to discuss is the kite integral family shown in figure 1.One of its subsectors is the equal-mass sunrise graph, whose ϵ 0 -part is known to involve complete elliptic integrals.In ref. [25], it was shown that there exists a choice of master integrals which fulfill differential equations in ϵ-form.The goal of this section is to extract the information about the elliptic functions from the differential equations matrix A(x, ϵ) and then reproduce the ϵ-form using our algorithm.

Definitions and differential equations
We define the kite integral family in D space-time dimensions and the equal-mass case as with In the following, we set m 2 ≡ 1 for simplicity as the dependence on this mass scale can be recovered from dimensional analysis.The integrals of the kite family then only depend on the kinematic variable Further, we introduce the dimensional regulator ϵ by setting D = D 0 − 2ϵ with D 0 being an even integer.In particular, we will be interested in the case D 0 = 4.
Using FIRE6 [46] and LiteRed [47], we find that there are a total of eight master integrals, which we choose as Note that the normalization of f 8 was chosen s.t.its integrand has a dlog-representation with constant leading singularity in D = 4 [12][13][14].Integrals of this type are expected to evaluate to pure functions [7,13], which makes f 8 a suitable initial integral for our algorithm.

Ansatz for the ϵ-form
To construct an ansatz according to (2.8), we first analyze the differential equations for the basis given in (3.4).We find that they have poles in x at 0, 1, 9 and ∞.Therefore we also restrict the integration kernels a k 1 ,...,km (x) dx to have poles at these points only.
To get information on the class of functions that can appear in ϵ-form, we follow the procedure described in section 2.2 and solve the differential equations at ϵ = 0 in each of the sectors.We find that only the sunrise sector formed by f 5 and f 6 requires functions beyond rational and polylogarithmic functions.In particular, we find where the ellipsis indicate terms from subsectors.Note that, here we follow [25] and perform this analysis in D = 2 instead of D = 4 space-time dimensions.Since the integrals in different space-time dimensions can be related through dimensional recurrence relations [48,49], this will not change the information we are trying to extract, namely the required class of functions for the ansatz.Eq. (3.5) can equivalently be written as the following second-order differential equation: where Ψ 1,2 (x) indicate the homogeneous solutions for the first line in (3.5).A standard choice for the solutions is given in appendix A. In summary, solving eq.(3.5) requires the introduction of a function which satisfies eq.(3.6) and therefore we need to include it when constructing the ansatz for the ϵ-form also in D = 4 − 2ϵ dimensions.Further following section 2.2, we take our ansatz to consist of the integration kernels described in eq.(2.16), subject to the condition that they are linearly independent and free of double poles.For the latter condition it is important to know the behavior of Ψ 1 near all singular points, see appendix B. In addition, from looking at known ϵ-forms of univariate elliptic Feynman integral families in the literature (see also the other examples in this paper), we expect that a maximum degree of two in Ψ 1 should be sufficient for the kite integral family.

The resulting ϵ-form
Using the initial integral f 8 = g 8 together with the ansatz discussed in the last section is sufficient input for our algorithm to reproduce the ϵ-form first found in [25]: where the integration kernels a i,j are given by the set Note that this is only a subset of all kernels used in the ansatz, meaning that some of the constant coefficient matrices m i,j were determined to be vanishing by our algorithm.As noted before, the transformation to the basis ⃗ g can be downloaded together with our public implementation, see section 4.

Two-loop non-planar triangle integral with internal masses
As a second application, we consider the integral family defined by the non-planar twoloop three-point graph depicted in figure 2. It has been analyzed in [50] by means of the differential equations method where it was shown that the top topology gives rise to two master integrals that cannot be written in terms of multiple polylogarithms.Instead, the expressions presented for their finite pieces involved integrals over products of complete elliptic integrals of the first kind and polylogarithmic functions.We note that an ϵ-form for the two integrals of the top sector was found in appendix A of [41].Here, we wish to demonstrate that our algorithm, supplied with a suitable ansatz, yields an ϵ-form for the full differential equations (including subsectors) that leads to similar iterated integrals order-by-order in ϵ as were found for the kite family.The propagators for this example are (3.9) For convenience, we set m 2 ≡ −1 from now on.The external momenta satisfy p 2 1 = p 2 2 = 0 and we define the kinematic variable (3.10) The graph from figure 2 corresponds to the sector G 1,1,1,1,1,1,0 .All integrals in this family can be expressed in terms of eleven master integrals, two of which lie in the top sector.
Reformulating the homogeneous part of the top-sector differential equations in D = 4 as a second-order differential equation for the scalar integral, we obtain Consequently, as for the kite family, we have to include one of the two solutions Ψ 1,2 (x) in our ansatz.For simplicity, we take the solution Ψ 1 (x).
Analyzing the diagonal blocks of the differential equations for the nine subsector master integrals, one finds that their solutions contain the two square-roots Hence, as was argued in section 2.2, our ansatz has to be rational not only in Ψ 1 (x), but also r 1 and r 2 . 3sing the same constraints as for the kite family in section 3.1.2,we can now write down a finite amount of terms for the ansatz.As the initial integral, we expect that the scalar integral divided by Ψ 1 , i.e.
is a good choice.This is based on the observation that in D = 4 the maximal cut of this integral evaluates to a pure function when integrated on a particular contour [41,51,52].
Applying the algorithm discussed in section 2.1, one encounters a particularity of this integral family: The sector G 1,1,0,1,1,1,0 admits one master integral that decouples completely from the remaining system of differential equations, such that they split into two separate problems.In principle, this would not pose a problem for our algorithm as we could bring them to ϵ-form simultaneously by taking a second initial integral from this sector.However, in this case, the differential equation for G 1,1,0,1,1,1,0 is already trivial.Therefore, we just discard it and proceed with the remaining ten master integrals.An ϵ-form is then easily obtained, where the independent integration kernels are given by the set This ϵ-form completes the result for the nine subsector integrals in [50] and the canonical differential equations for the maximal cut of the top sector presented in appendix A of [41], by also providing a transformation for the off-diagonal blocks.Further, this example demonstrates nicely how square-roots are handled easily by the algorithm.

N3LO Higgs-production phase-space integrals
As a third example, we consider some elliptic phase space integrals which appear in the calculation of the partonic coefficient functions for Higgs production at N 3 LO in QCD [53].The associated graph is depicted in figure 3 and the propagators are where we use the notation p i 1 ...in = p i 1 +. ..+p in .Note that D 1 , . . ., D 4 are cut propagators, which means that integrals with non-positive a 1 , . . ., a 4 are zero.The kinematics is s = (p 1 + p 2 ) 2 , p 2 1 = p 2 2 = 0 and we set x = m 2 h and s = 1.We work in D = 4−2ϵ space-time dimensions and consider the sector G 1,1,1,1,1,1,1,1,0,0,0,0 together with its four non-trivial subsectors.They give rise to a total of 19 master integrals.Examining the diagonal blocks of the differential equations at ϵ = 0, one finds in the top sector the second-order differential equation In [53] it was shown that the two independent solutions can be expressed in terms of complete elliptic integrals of the first kind.Following the same principles as before, an ansatz for the ϵ-form including one of the two, for concreteness Ψ 1 , is then written down easily and similarly to the previous section, we expect that a suitable initial integral is given by However, one finds that the first 19 derivatives of this integral are not linearly independent.Therefore, to construct an invertible Ψ-matrix, we need to supplement the derivatives of additional initial integrals.Through an integrand analysis or simply through trial and error, we find that already ten of the 19 basis integrals given by the IBP-reduction code FIRE6 are uniform weight integrals and therefore this is an easy task.For more information on the initial basis ⃗ f , we refer to the examples that come with our implementation, see section 4.
The independent integration kernels in the resulting ϵ-form are x (x 2 − 11x − 1) .
(3.18)In this section, we discuss a family of integrals that appears in the computation of the gravitational potential of two non-spinning binaries, see e.g.[54][55][56][57].The propagators in this case are

Three-loop gravitational potential integrals
where D 1 , D 2 and D 3 are again cut propagators.The kinematics is u Further, one sets γ = (x + 1/x)/2 to rationalize appearing squareroots and q 2 = −1 since this is the only dimensionful scale and therefore the dependence of the integrals on it is trivial.
We are considering the sector G 1,1,1,0,0,0,0,0,1,1,1,0,0,1,1 depicted in figure 4, which has three master integrals.Together with its subsector, there are a total of four master integrals.Note that the differential equations for the full family, which are needed for the computation of the gravitational potential, have around 60 master integrals and have likewise been brought into ϵ-form in [56] by complementing the algorithm presented in this paper with other methods.
Inspecting the differential equations for the scalar integral in D = 4, we find for the first time in this paper a third-order, rather than a second-order, Picard-Fuchs equation: However, one can easily verify (see e.g.[40]) that the solutions to this equation can be written as the products where Ψ 1,2 are the solutions of the following second-order differential equation: Similar to before, we take the scalar integral divided by Ψ1 = xΨ 2 1 as initial integral: G 1,1,1,0,0,0,0,0,1,1,1,0,0,1,1 . (3.23) However, because the Picard-Fuchs equation is of third-order, we find that the maximal power of Ψ 1 in the ansatz now has to be four instead of two. 4 The integration kernels in the resulting ϵ-form are We see that our algorithm is also applicable to the case of integrals that include functions satisfying higher-order Picard-Fuchs equations, as long as we can still make a reasonable guess for the ansatz.

Three-loop equal-mass banana graph
As the last example, we consider the three-loop equal-mass banana graph depicted in figure 5.This integral family has received a lot of attention due to being the natural extension of the sunrise graph to the next loop order.An ϵ-form has recently been achieved in [27], see also [40,58,59].The propagators are and we take the number of dimensions to be D = 2 − 2ϵ.The kinematics is p 2 = −(x − 1)(x−9)/x and m 2 = 1.For the banana graph, we consider the three master integrals of the sector G 1,1,1,1,0,0,0,0,0 together with the single master integral from its subsector.Therefore, there are a total of four master integrals.Setting ϵ = 0, we find the third-order differential equation for the scalar integral, similarly to the example in the previous section.Again, one can verify that the solutions can be written as where now Ψ 1,2 are the very same functions appearing in the sunrise graph in eq.(3.6).We therefore take as initial integral and construct our ansatz accordingly.However, given our ansatz and the initial integral f 1 , we find that there is no solution to (2.10), suggesting that either (3.28) is not a suitable integral or that our ansatz needs to be supplemented with additional functions.Indeed, the known ϵ-form in [27] involves new integration kernels which we do not find through the procedure described in sections 2.2 and 3.1.2.In particular, they involve a new function, called F 2 in [27], which satisfies the differential equation5 techniques necessary for their calculation to a similar level as for polylogarithmic Feynman integrals.In this paper, our aim was to contribute to this by extending the algorithm of [23] to include functions satisfying higher-order Picard-Fuchs equations.In particular, we described how to extract information about these functions and proposed a suitable way of making an ansatz for the precise form of the canonical differential equations.This ansatz has two main inputs: Firstly, it is motivated by known ϵ-forms in the literature, and by the expected target class of iterated integrals.Secondly, as in the polylogarithmic case, it uses integrand analysis (considering multiple residues, or generalized cuts) to find a suitable 'initial' integral, whose Picard-Fuchs equation is crucial input for our algorithm.The latter then finds a complete basis of canonical integrals, assuming the ansatz for the canonical differential equations contained all relevant terms.We discussed several state-ofthe-art examples, including new cases, such as the full form of the ϵ-factorized differential equations for the massive form factor integrals in section 3.2, and the Higgs phase space integrals in section 3.3.Our work opens up several possible directions for further investigation: Multivariate kinematics: Although we focus on examples with only a single variable x, the algorithm applies to multivariate cases as well, see [23] section 3.4.A natural first example of this is the unequal-mass sunrise integral family, for which an ϵ-form has been computed in [26].The integration kernels then also involve functions satisfying an inhomogeneous second-order differential equation that can be handled by our algorithm in a similar way as F 2 of section 3.5.
Two elliptic curves: There are known examples of Feynman integral families where two of the diagonal blocks give rise to two different kinds of elliptic functions [28,61].The offdiagonal blocks can then simultaneously depend on both of these sectors, which complicates the process of finding an ϵ-form.We are confident that our algorithm should be able to reproduce the known results [28].However, it would certainly be interesting to motivate an ansatz without any information on the known ϵ-form, because this should then also allow us to apply the algorithm to other problems of similar nature.
Higher-order Picard-Fuchs equations: We have already seen in sections 3.4 and 3.5 that differential equations involving a third-order differential equation at ϵ = 0 can be treated in a similar way as for second-order differential equations.However, both of these examples eventually turned out to involve elliptic functions only.It would be interesting to see how our algorithm can be applied to examples that involve functions that go beyond elliptic integrals (see e.g.[62]).
Modularity and numerical integration: For some of our examples, the new integration kernels involve square-roots, or their denominators have algebraic roots.Therefore, it would be interesting to study their behavior under modular transformations, with the idea of relating them to functions better suited for numerical integration, see e.g.[32].
The series expansions of F 2 can also be derived either by using the method of Frobenius or Wasow in combination with variation of constants, or by starting from the explicit integral representation.The result depends on which representatives are chosen for the two periods If we always choose ψ 1 = ψ x 0 ,1 for every singular point x = x 0 , the expansions are simple power series: Other choices lead to very complicated expressions involving logarithms.For example

Figure 1 .
Figure 1.The kite integral family.Thick lines denote massive propagators.If one pinches the two massless lines, one recovers the sunrise graph which is the source of elliptic integrals.

Figure 2 .
Figure 2. A non-planar two-loop family giving rise to elliptic integrals.The thick lines denote massive propagators whose masses are taken to be equal.

Figure 3 .
Figure 3. Elliptic sector in the phase-space integrals for Higgs boson production at N 3 LO in QCD.The thick line denotes the massive Higgs line.Lines crossing the dashed line denote cut propagators.

Figure 4 .
Figure 4.A three-loop integral sector that gives rise to elliptic functions in the gravitational potential of non-spinning binaries.Double lines denote linear propagators.

Figure 5 .
Figure 5.The three-loop equal-mass banana graph.The mass of the internal lines is denoted by m.