All orders structure and efficient computation of linearly reducible elliptic Feynman integrals

We define linearly reducible elliptic Feynman integrals, and we show that they can be algorithmically solved up to arbitrary order of the dimensional regulator in terms of a 1-dimensional integral over a polylogarithmic integrand, which we call the inner polylogarithmic part (IPP). The solution is obtained by direct integration of the Feynman parametric representation. When the IPP depends on one elliptic curve (and no other algebraic functions), this class of Feynman integrals can be algorithmically solved in terms of elliptic multiple polylogarithms (eMPLs) by using integration by parts identities. We then elaborate on the differential equations method. Specifically, we show that the IPP can be mapped to a generalized integral topology satisfying a set of differential equations in ϵ-form. In the examples we consider the canonical differential equations can be directly solved in terms of eMPLs up to arbitrary order of the dimensional regulator. The remaining 1-dimensional integral may be performed to express such integrals completely in terms of eMPLs. We apply these methods to solve two- and three-points integrals in terms of eMPLs. We analytically continue these integrals to the physical region by using their 1-dimensional integral representation.


JHEP01(2019)169 1 Introduction
The evaluation of Feynman diagrams is a crucial ingredient of particle physics calculations. Analytic evaluation is required up to high precision in order to make predictions for processes at the LHC, or other particle colliders. The evaluation of Feynman diagrams becomes increasingly complicated as the number of loops and scales increases. Many techniques have been developed to deal with this complexity. Modern methods focus on scalar Feynman diagrams, which span the diagrams for more general (gauge) theories. This result follows from applying Passarino-Veltman reduction [1] and integration-by-parts (IBP) identities [2][3][4][5] to rewrite non-scalar diagrams into scalar ones. In the remainder of the paper we will always let the term Feynman diagram refer to scalar Feynman diagrams, which we also call Feynman integrals.
Many Feynman integrals can be analytically computed in terms of special functions known as multiple polylogarithms (MPL) [6,7]. This could be considered the "simplest" class of functions that one may hope to consider, as for example the massless on-shell (1-loop) bubble, triangle, box and pentagon integrals are expressible in terms of MPLs up to all orders in the dimensional regulator ǫ. Upon considering more complicated integrals, by adding external massive legs, internal massive lines, or considering higher loop diagrams, one quickly encounters integrals that are not expressible in terms of MPLs, and more general classes of functions need to be considered.
From the traditional viewpoint Feynman integrals are defined as momentum space integrals: one integrates over a 4-dimensional space for each internal momentum. Although many such integrals are divergent, one may adopt a dimensional regularization scheme such that the dimension is taken to be d = 4−2ǫ, where ǫ is the dimensional regulator. Different integer dimensions than 4 may also be considered in this scheme. Unfortunately from the momentum-space viewpoint Feynman integrals are difficult to calculate both analytically and numerically. To get a better handle on them it is necessary to study things from a different viewpoint. Arguably the two most successful but different starting points are the differential equations method [8][9][10][11][12] and the direct integration of the scalar parametrization of Feynman integrals [13,14]. We give a short discussion of these two methods next.
The differential equations method has received significant attention in the last decade. In this method collections of Feynman integrals are considered which form a basis for the system of IBP relations that is associated with their topology. 1 Such integrals are referred to as master integrals. The set of master integrals may be minimized by using symmetry relations. For a basis of master integrals one may define a closed system of differential equations under differentiation with respect to external scales, like the Mandelstam variables and the masses. The study of such a system has been systematized by the introduction of the so-called canonical basis of master integrals [15], which is conjectured to exist for polylogarithmic integrals. For such a basis of master integrals, the differential equation matrix is ǫ-factorized, where ǫ is the dimensional regulator, and the differential equation matrix is in dlog-form. The general solution of such a system can be written JHEP01(2019)169 as a path-ordered exponential, which gives the individual terms in the ǫ-expansion as iterated integrals. Furthermore, if the arguments of the dlog's are rational, the resulting Chen-iterated integrals can be directly expressed in terms of multiple polylogarithms. The problem of finding the canonical basis -when it exists -has been automated in a number of software packages [16][17][18]. Another recurring challenge in the differential equations method is to perform the IBP reduction to obtain a closed form for the differential equations. Specialized software packages are available to perform these reductions [19][20][21], but all run into computational limits for sufficiently complicated topologies.
The approach of directly integrating the scalar parametrization of Feynman integrals has seen major advancements by the works of (among others) Brown and Panzer, see for example [13,14]. In particular the Maple package HyperInt [22] automates the integration of multiple polylogarithmic functions. The direct integration method works as follows. First one computes the parametric (Feynman) representation of a Feynman integral, which will be discussed in section 2. Starting from this representation one attempts to integrate one integration parameter at a time. In the ideal scenario there exists an integration sequence such that all integrations can be done in terms of multiple polylogarithms. In that case the Feynman integral is called "linearly reducible", as it requires a linear factorization of the set of polynomials occurring in each integration step. If instead a linear factorization can only be performed at the expense of introducing algebraic roots that contain remaining integration variables, one will generally not be able to perform these integrations in terms of MPLs. (Although in some cases a change of variables may be obtained to transform away the square roots.) The property of linear reducibility can be investigated without explicitly performing the integration, by considering if the compatibility graph associated to the parametric representation is linearly reducible [13,14].
The next challenge in the computation of Feynman integrals is to get a good grasp on the structure of integrals that can't be expressed in terms of MPLs. It seems that the next non-trivial class are the so-called elliptic Feynman integrals, which have seen a lot of interest in the last years . Elliptic integrals have traditionally been (loosely) identified as integrals for which the maximal cut, or the maximal cut of one of their subtopologies, is an elliptic integral.
Considerable progress has been made in the last years in the understanding of elliptic Feynman integrals. We discuss a few different starting points for solving elliptic Feynman integrals next. In [38] a non-planar triangle topology is considered for which the top sector is elliptic. In that case one may first find and solve a canonical dlog basis for the subtopologies, in terms of multiple polylogarithms. The solution for the elliptic Feynman integrals in the top sector may then be written down by the method of variation of parameters. This requires solving the homogeneous equations of the integrals in the top-sector. The problem of finding these homogeneous solutions can in turn be tackled by computing the maximal cuts of the integrals in the top sector. Since these maximally cut integrals have vanishing subtopologies, they solve the homogeneous part of the differential equation for the uncut integral (the inhomogeneous terms in the differential equation are subtopologies). The computation of maximally cut Feynman integrals has furthermore been aided by recent explorations of the Baikov parametrization [48][49][50].

JHEP01(2019)169
The equal-mass sunrise, which has tadpole integrals as subtopologies, has also been extensively studied in this way, see for example [25,28,30]. In more recent work it has been shown that an ǫ-factorized form for the differential equations of the kite integral family, which contains the equal mass sunrise, may be obtained by allowing for non-algebraic integration kernels in the differential equations matrix. The resulting solution is then given in terms of iterated integrals of modular forms [46]. A planar double box integral relevant to top-pair production with a closed top loop, involving several elliptic sub-sectors, has furthermore been computed by considering a differential equation linear in ǫ [47,51].
On the more formal side significant progress has been made in the understanding of iterated integrals over an elliptic curve. More specifically, all iterated integrals over an elliptic curve can be expressed in terms of special functions called elliptic multiple polylogarithms (eMPLs) [44,52,53]. Moreover, their functional identities can be systematically studied by using a suitable elliptic generalization of the symbol map [54]. While this class of functions seems to be the natural candidate to express a large class of elliptic Feynman integrals, it is still unclear how and when this representation can be obtained.
In the first part of the paper we investigate dimensionally regulated elliptic Feynman integrals by direct integration of the Feynman parametric representation. In so doing we identify a class of Feynman integrals that we call linearly reducible elliptic Feynman integrals. These integrals can be algorithmically solved up to arbitrary order of the dimensional regulator in terms of 1-dimensional integrals over a polylogarithmic expression, the so-called inner polylogarithmic part (IPP). More specifically, the direct integration approach requires that the polynomials in each integration step factor linearly with respect to the next integration variable (Feynman parameter). However, in the elliptic case, the factorization introduces irreducible square roots at some integration step and no further integration can be done thereafter within the class of multiple polylogarithms. Nevertheless we have observed empirically that in many cases the irreducible square roots appear only when the integration with respect to the last parameter has to be performed. We call this class of elliptic Feynman integrals linearly reducible. When the IPP depends on one elliptic curve and no other algebraic functions this class of Feynman integrals can be algorithmically solved in terms of eMPLs by using, e.g., integration by parts identities [44].
While in many cases the direct integration approach is convenient as, e.g., no boundary conditions need to be computed and results at relatively low orders of the dimensional regulator can be easily obtained, it might be impractical to apply it when higher orders need to be considered. One reason is that the size of the expressions usually grows very rapidly and the algebraic manipulations involved become cumbersome. Moreover, the analytic properties of the answer to all orders are not immediately manifest in this approach. For these reasons we believe that a better understanding of the differential equations method applied to Feynman integrals beyond multiple polylogarithms is highly desirable. In the second part of the paper we show that the IPP of linearly reducible elliptic Feynman integrals can be mapped to a generalized integral topology satisfying a set of differential equations in ǫ-form. The mapping is obtained by applying the Feynman trick to a suitably chosen pair of propagators in the momentum space representation of the Feynman integral under consideration. Remarkably, for all the examples we considered the canonical JHEP01(2019)169 differential equations matrix can be directly expressed in terms of linear combinations of eMPL integration kernels, and the solution of these equations in terms of eMPLs becomes elementary. Once the IPP is expressed in terms of eMPLs the remaining 1-dimensional integral can be easily performed to express linearly reducible elliptic Feynman integrals completely in terms of eMPLs.
The remainder of the paper is organized as follows. In section 2 we review the properties of the parametric Feynman representation of Feynman diagrams. In particular we discuss the Cheng-Wu theorem, which is needed to find a linearly reducible integration order, and we review the direct integration algorithm and some of its properties. In section 3 we discuss linearly reducible elliptic Feynman integrals and some properties of the class of 1-fold integrals that arises in their computation. In section 4 we explain how to set up differential equations for the inner polylogarithmic part of a linearly reducible elliptic Feynman integral, which involves an application of the Feynman trick and is thereafter analogous to the polylogarithmic case. We then move on to section 6, where we treat a few examples of linearly reducible elliptic Feynman integrals. In section 6.1 we discuss the unequal mass sunrise in depth. We study it by direct integration, and as a solution from the system of differential equations for its IPP. We also derive the analytic continuation of the first master integral to the physical region. A "triangle with bubble" integral is discussed in section 6.2 for which we perform a similar analysis. In section 7 we solve an elliptic non-planar triangle integral relevant for Higgs+jet production. This integral is linearly reducible, however the IPP depends on multiple algebraic curves and further development of the methods studied in this paper will be required. We end our discussion in section 8, where we reflect on the results obtained and give an outlook for the future.

Parametric representation of Feynman integrals
The defining momentum space integral representation of scalar Feynman diagrams admits the form: where {s} = {p} ∪ {m} schematically denotes the dependence on the set of external momenta and masses, which is left implicit on the right hand side of the equation. Furthermore, N is a normalization constant picked by convention, l denotes the number of loops, n denotes the number of propagators and d is the dimension of the Minkowskian space-time.
The convention D i = −q 2 i + m 2 i is used, where q is the momentum flowing on the i-th propagator and m is the mass on the propagator. We assume that the exponents a i of the propagators are positive integers.
A scalar parametrization can be found as follows. First one uses the Schwinger trick to rewrite every propagator as: 2 (2.2)

JHEP01(2019)169
The momentum integrals of eq. (2.1) become Gaussian integrals which can be performed and result in the so-called Schwinger parametrization. From there, one may perform a change of variables α i → ηα i , under the constraint n i=1 α i = 1. The parameter η can be integrated out, which yields the Feynman parametrization: where U and F are called Symanzik polynomials, which are related to determinants taken during the Gaussian integration, but have an alternative interpretation from considering the set of spanning trees T (G) and spanning 2-forests 3 F 2 (G) of the Feynman diagram G: where s (T 1 ,T 2 ) is defined as the square of the external momentum travelling in between the components of the 2-forest. For a review of these concepts see for example [55,56]. The region ∆ is defined as The scalar parametrization of eq. (2.3) is called the Feynman parametrization and the integration parameters are called Feynman parameters. We will also refer to it as the parametric representation. In the remaining sections we generally choose the normalization constant (2.5) to remove the prefactor in the Feynman representation. We will often consider the external kinematics in the Euclidean region, which is determined by the condition that F ≥ 0 on the whole domain of integration.

The Cheng-Wu theorem
The Cheng-Wu theorem specifies various ways that one might perform the integration over the Feynman parametrization. Importantly, different applications of the theorem may be helpful in finding a linearly reducible integration order. The Cheng-Wu theorem states that a projective (Feynman) integral over ∆ has the same value when integrated over the domain where S ⊆ [1, n] is a nonempty set of integers [57]. Projectivity refers in this context to the property that the integrand should remain invariant under the rescaling: 3 A 2-forest is a disjoint union of 2 trees.

JHEP01(2019)169
If one starts with an integral over ∆ which does not remain invariant under eq. (2.7), the following projective transform can be performed to obtain such an integral which has Jacobian ( n j=1 α ′ j ) −n [58]. Note that this change of variables keeps the integration domain ∆ invariant and is in fact only a proper change of variables because we integrate over ∆, as it manifestly sets the sum of the integration parameters to 1. As an illustrative example consider the beta function: Eq. (2.9) starts with the usual definition of the beta function, upgrades it to a 2-dimensional integral over ∆, projectivizes the integrand, and lastly uses Cheng-Wu to achieve a different representation of the integral. (Of course in this trivial example the same result would have easily been obtained with the Möbius transformation α 1 → α ′ 1 /(α ′ 1 + 1).) To integrate the scalar parametrization one seeks a linearly reducible integration order, which is discussed in the upcoming section. As a rule of thumb it is often sufficient to apply Cheng-Wu before performing any integrations, and to let one of the Feynman parameters, say α i , go to 1 by picking S = {i}. The other parameters are then integrated from 0 to infinity. It is sometimes convenient to factor out the Cheng-Wu parameter: (2.10) One can then focus on integrating out Feynman parameters successively in the bracketed part. Projectivity of the integrand is usually manifestly preserved after performing these integrations. Suppose one ends up with the following expression after performing a number of integrations: where there are k < n − 1 non-trivial Feynman parameters remaining, labeled by a set S ′ , and where the integrand is projective. By the Cheng-Wu theorem we can turn this into: where S ′′ ⊆ S ′ ∪ {i}.

Direct integration, linear reducibility and all orders statements
First we remind the reader of the definition of (Goncharov) multiple polylogarithms (MPLs). They are the following recursively defined functions: G( a n ; z) = G(a 1 , . . . , a n ; z) = z 0 dt t − a 1 G(a 2 , . . . , a n ; t) , where a 1 , . . . , a n , z are complex variables. The recursion is ended at n = 0 where one lets by convention: As a general feature of iterated integrals, MPLs obey the shuffle product: where the set of shuffles of a n and b m denoted a n ⊔⊔ b m may be understood to contain all permutations of the sequence ( a, b) that preserve the ordering of the individual vectors. A (small) technical complication in the definition of multiple polylogarithms is that a divergence at the basepoint 0 occurs when a n = 0. One may adopt the definition: 16) to deal with the divergent case where all n parameters a i are equal to zero. Cases with a n = 0 and at least one a i = 0 can then be dealt with in a consistent manner by rearranging parameters using the shuffle product, and using eq. (2.16). A pedagogical review of multiple polylogarithms and their functional identities is given in [59]. The parametric representation discussed in the previous sections can be integrated, e.g., by using the computer program HyperInt [22]. We sketch a few ideas underlying the algorithm next, in order to illustrate the concept of linear reducibility. First one performs a series expansion of the integrand of eq. (2.3) on the dimensional regulator. Assuming that the Feynman integral is finite in the integer dimensiond, we let d =d − 2ǫ and find: where the coefficients are: It is clear that for even dimensionsd and integer powers of the propagators the resulting integrand is a polylogarithmic expression without algebraic terms. The remaining integrations to be performed, at each order in ǫ, take the schematic form:

JHEP01(2019)169
where k ∈ {1, . . . , n}. We aim to perform these integrations in such a way that at each integration step the integrand f k−1 (α k , . . . , α n ) is a polylogarithmic expression. More precisely, we require that the integrand at each integration step is a combination of multiple polylogarithms with prefactors and arguments that are rational functions of the remaining integration parameters. Now suppose that f k−1 (α k , . . . , α n ) is a polylogarithmic expression. Then it depends on a set of irreducible polynomials in the remaining integration parameters which we denote by P (k−1) (α k , . . . , α n ). A requirement for the integral in eq. (2.19) to be polylogarithmic again is that all polynomials in P (k−1) (α k , . . . , α n ) are linear in α k . If at each step of the integration the set of polynomials P (k−1) (α k , . . . , α n ) is linear in α k , we have found that α 1 , . . . , α n is a so-called linearly reducible integration order.
Each integration can then be performed in terms of multiple polylogarithms along the following lines: 1. express f k−1 (α k , . . . , α n ) as a combination of multiple polylogarithms of argument α k , 2. find a primitive F k (α k , . . . , α n ) such that ∂ α k F k (α k , . . . , α n ) = f k−1 (α k , . . . , α n ), To search for a linearly reducible integration order one can enumerate over all possible integration sequences. Luckily the set of polynomials P (k) (α k+1 , . . . , α n ) at each integration step can be exposed from a so-called compatibility graph without performing the actual integration, as introduced in [13]. While we refer the reader to that reference for further details, an important consequence of the compatibility graph method is that the polynomials P (k) (α k+1 , . . . , α n ), k ∈ {1, . . . , n} are independent of the order in ǫ we are considering. It should however be noted that at leading order in ǫ the exponent of one of the Symanzik polynomials may become 0 for special configurations of the dimension and powers of the propagators, and hence one may find an integration sequence that does not work at higher orders in ǫ. In section 6 we solve a number of finite linearly reducible elliptic Feynman integrals up to and including the order ǫ 1 , and it is understood that this yields linear reducibility up to all orders. In some cases one finds that applying the Cheng-Wu theorem with one Feynman parameter set to 1 does not lead to a linearly reducible integration order. A nontrivial application of the Cheng-Wu theorem may then occasionally lead to a linearly reducible integration order. This is the case for the non-planar triangle integral in section 7. In other cases the situation may be worse, and a change of variables is needed at some point during the integration. Some more discussion on this topic, and explicit examples of changes of variables are discussed in [14,60].
The previous story applies when the Feynman integral is expressible in terms of multiple polylogarithms. To our knowledge there is no known example of an elliptic Feynman integral which can be expressed in terms of MPLs, and it is believed that such a representation is not possible for elliptic Feynman integrals in general. That means in particular that no linearly reducible integration order exists for these integrals. We therefore define JHEP01(2019)169 linearly reducible elliptic Feynman integrals as elliptic Feynman integrals that are linearly reducible if one excludes the last integration. In particular, we allow the next-to-last integration to be performed at the expense of introducing algebraic terms of the last integration parameter. The final expression of a linearly reducible elliptic Feynman integral reads, at a generic ǫ-order: where f n−1 (α n ) is a polylogarithmic expression with algebraic coefficients and arguments, depending on the elliptic curves of the problem. In the framework of the direct parametric integration, linearly reducible elliptic Feynman integrals are the simplest instance of elliptic Feynman integrals.

Structure of linearly reducible elliptic Feynman integrals
In the previous section linearly reducible elliptic Feynman integrals have been introduced which are, to all orders in ǫ, expressible as 1-fold integrals. These 1-fold integrals take the following schematic form: where the sum over i denotes a generic collection of terms grouped by factors A i which are algebraic functions in x, and MPL i (x) denotes a polylogarithmic term with algebraic arguments. Here we investigate some properties of these 1-fold integral, and sketch some general strategies that may be employed to write eq. (3.1) in terms of a minimal class of integrals. From the direct integration point of view the algebraic dependence in A i arises from forcing a linear factorization of the polynomials in the previous integration step. For the upcoming examples of the unequal mass sunrise, and the triangle with bubble, we find that the only algebraic dependence of the inner polylogarithmic part is on 1 elliptic curve. In that case one may write: such that y(x) 2 = P (x) defines an elliptic curve, i.e. P (x) is an irreducible cubic or quartic polynomial, and where R i is a rational function in its arguments. One may furthermore factorize the dependence on the elliptic curve such that one has: where R i,1 (x) and R i,2 (x) are rational functions in x. This can be achieved in the following manner. Firstly one may absorb any even power of y(x) in the rational part. Furthermore, for denominators of the type 1/(S 1 (x) + S 2 (x)y(x)) k , where S 1 (x) and S 2 (x) are some polynomials in x and k is some positive integer, one may multiply by the conjugate

JHEP01(2019)169
observe that the new denominator is a polynomial, and expand out the numerator, again absorbing all even powers of y(x) in the rational part in x. Lastly one may use the relation y(x) = P (x)/y(x) to obtain a representation of the form of eq. (3.3). One may furthermore partial fraction a rational term R(x) such that: is a polynomial in x, and N i,j and M j are complex coefficients that do not depend on x. From now on we will shorten the notation y(x) to y. By the previous arguments we may reduce the integrand to the following cases: where k is a positive integer, β is a constant with respect to x, and MPL(x, y) is a polylogarithmic expression. We will refer to the algebraic factors, without the polylogarithmic term MPL(x, y), as integration kernels. Note that splitting up the integral in eq. (3.1) in terms of 1-fold integrals of the type of eq. (3.6) requires that the individual contributions are finite. We will ignore the issue of regularization of the individual contributions for now, and provide results in terms of a single one fold integral that contains all contributions. Furthermore, to save space we will use the shorthand subscript notation: G(a 1 , . . . , a n ; 1) = G a 1 ,...,an , (3.7) later on in the text to present the polylogarithmic terms. In the upcoming examples we will find that we are able to pick a basis of master integrals such that we only encounter kernels with k = 1 and with an elliptic curve that is quartic. Nonetheless, in general one might expect other integration kernels to show up as well, and we show next that it is possible to reduce kernels of the type with k ¿ 1 to the case with k = 1 by employing IBP relations. This is similar to the treatment in [44], where these kinds of IBP identities are considered to provide an integration algorithm for elliptic polylogarithms multiplied by rational functions. In our treatment we keep the factor multiplying the algebraic integration kernel explicitly polylogarithmic, and we will work with a quartic elliptic curve whose roots will be denoted by a 1 , . . . , a 4 . Besides reducing the integration kernels to the case with k = 1, it is possible to relate kernels of the form dx/((x − a i )y). There are 4 such kernels, and one may trade them using IBP identities for one of the following:

JHEP01(2019)169
where we note that the polylogarithmic part MPL(x, y) that multiplies the kernels is affected by the IBP relations. The final result that is obtained after reducing the set of integration kernels is not necessarily unique. For example, the inner polylogarithmic part is still subject to the usual functional identities between multiple polylogarithms. Furthermore, a kernel of the type dx/(x − β) may be exchanged for different ones using the IBP identity: where the new kernels depend on the precise form of MPL ′ (x, y).
Next we provide the explicit IBP relations that may be used to reduce the set of integration kernels to the cases with k = 1. First we consider the following relation for k > 1: where we note that: from which it is clear that the power of x in the numerator is reduced for each term in the indefinite integral on the right hand side that carries a factor MPL(x). Note that there is a term carrying a factor x k+1 , but this is multiplied by the derivative of MPL(x), which has its weight reduced by 1. This allows for an inductive procedure that terminates at weight 0. Similarly we may derive the following relation for k > 1: (3.12) We remark that partial fractioning a term of the form ( , and a piece that carries a factor 1/(x − a i ), see eq. (3.5). Hence we may safely iterate eq. (3.12) to reduce the power of (x − c) in the denominator, up to polylogarithmic terms of lower weight. Lastly, we provide the following relation that may be used to trade the kernel dx/(y(x − a 1 )) for x 2 dx/y, up to polylogarithmic terms of lower weight:

JHEP01(2019)169
where we used the notation a ij = a i − a j . One may obtain similar relations for the kernels dx/(y(x − a j )), j = 2, 3, 4, by cyclically permuting the labels of the roots: a i → a i+1 . This way one may remove every kernel of the type dx/(y(x − a i )), at the expense of introducing a kernel x 2 dx/y. One may rearrange eq. (3.13) afterwards to obtain an expression that contains just the kernel dx/(y(x − a 1 )). Lastly we have the relations: of which the right hand side in both cases involves a piece that has been integrated, and an indefinite integral that contains terms of lower weight.

Differential equations for the inner polylogarithmic part
In the previous sections we have considered linearly reducible elliptic Feynman integrals from the viewpoint of the direct integration method. In section 4.1 we will show that the inner polylogarithmic part of these integrals can be mapped to a (generalized) Feynman integral topology that arises from an application of the Feynman trick to two propagators. This topology can be studied in momentum space, where it is easy to derive IBP relations and setup a system of differential equations for its master integrals. We then review the differential equations method in section 4.2. In particular, we discuss the canonical basis of differential equations and how, in the upcoming examples, it can be used to algorithmically solve the IPP in terms of eMPLs. The full integral is then solved in terms of eMPLs by performing the remaining 1-fold integral over the IPP (see section 6). The resulting approach essentially bridges a gap between the direct integration method and the differential equations method. One can either find the IPP by direct integration of the Feynman parametrization, or alternatively one can find it by solving a canonical system of differential equations for its corresponding topology. However, in the direct integration approach the solution of a given integral in terms of eMPLs involves first integrating the IPP in terms of MPLs, and then iteratively writing the polylogarithms as an integral over their derivative. The complexity of this approach usually grows quickly with the order of the dimensional regulator. On the other hand in section 6 we show that by using the canonical differential equations method for the IPP it is possible, in some cases, to solve the full integrals up to arbitrary order of the dimensional regulator in a fully algebraic manner, since the relevant integration kernels coincide with the ones defining the eMPLs.

The Feynman trick
In the following treatment we consider the inner polylogarithmic part with respect to the last integration parameter α n−1 , which is a generic choice as we have the freedom to relabel variables. We start by considering a general topology I a 1 ,...,an . First write down

JHEP01(2019)169
the Feynman parametrization, and apply the Cheng-Wu theorem to put α n = 1: where we denote the inner polylogarithmic part with respect to the last integration on α n−1 as IPP (n−1) . The Feynman trick tells us that: Γ (a n−1 + a n ) Γ (a n−1 ) Γ (a n ) Inspired by this, we consider a new topology that contains a generalized propagator of the form α n−1 D n−1 + D n . In this topology α n−1 is an external scale, and to avoid confusion we'll denote the Feynman parameters of this topology with a tilde (α i ). We define: Like before a normalization factorÑ is included to remove an overall prefactor from the Feynman parametrization: From eq. (4.2) it is clear that this yields: ..,a n−2 ,a n−1 +an . (4.5) Next we show explicitly that: ..,a n−2 ,a n−1 +an .
Setting up the Feynman parametrization yields: I a 1 ,...,a n−2 ,a n−1 +an = where the Cheng-Wu theorem has been applied to put the Feynman parameterα n−1 (which corresponds to the generalized propagator) to 1. Next let U (α 1 , . . . , α n ) and

The differential equations method
Next we remind the reader of some points that are relevant for the differential equations method. Firstly one requires a reduction of the integrals in the topology in terms of a set of master integrals, which we denote by B = (B 1 , . . . , B k ). The master integrals are independent with respect to all IBP relations. (Such IBP relations are most easily generated from the momentum space picture of the integrals.) Furthermore, generally one also takes into account symmetry relations. A set of master integrals, and the reduction of the remaining integrals in the topology -up to some finite bound on the propagator exponents -may be obtained using programs such as LiteRed, FIRE and KIRA [19][20][21]. We will make use of the C++ version of FIRE5 which seemingly has no problem in dealing with combined propagators that are obtained from the Feynman trick. One may write down a closed form system of differential equations for B with respect to each external scale s j : whereÃ j is a matrix whose elements depend on the external scales and the dimension. For polylogarithmic topologies a basis can be found in a "canonical" d log ǫ-factorized form [15] where the matrices satisfy:Ã

JHEP01(2019)169
such that A has no more dependence on the dimension d =d − 2ǫ, whered is an integer and the dimensional regulator is ǫ. The set of "letters" A consists of rational or algebraic functions of the external scales. Lastly, A l is a matrix with integer coefficients. The differential equations for each s i may now be combined: Differential equations in canonical form have two important properties. Firstly, since ǫ is factored out one may write the general solution of the equation in terms of a path-ordered exponential: which order-by-order in ǫ expresses the result in terms of iterated integrals: where γ is a path with domain [0, 1] in the space of external invariants, and where we have a boundary term B boundary , which is B evaluated at the point in kinematic space given by γ(0). We furthermore denote the ǫ-expansion of the boundary term by: which we assume to be finite. Note that one may obtain a set of master integrals that is finite as ǫ → 0 by multiplying all the master integrals by a power of ǫ. If the letters are rational functions, and one is able to find a boundary term, eq. (4.14) allows us to directly write down the master integrals in terms of MPLs order by order in ǫ.
The second important result is that a canonical form differential equation provides the symbol of the master integrals, given that they are uniformly transcendental and that we have the leading coefficient of B in ǫ. In particular, one finds that: , and where R is an operator that reverses the ordering of the tensor product: In the upcoming sections we will also use the differential equations method to find results in terms of E 4 -functions [44,45], for examples that depend on a single quartic elliptic curve. This will be done by rescaling the last integration parameter, so that it runs from 0 to 1. We will denote the rescaled parameter as x ′ , and consider a system of differential equations with respect to x ′ :

JHEP01(2019)169
where B will be a canonical basis for the inner polylogarithmic part. The solution in terms of E 4 -functions will follow because the entries of the matrix ∂A/∂x ′ will turn out to correspond to linear combinations of integration kernels of E 4 -functions. The particular kernels that show up in the upcoming sections are presented here for completeness: For the definitions of the other integration kernels we refer to [44]. We note that c 4 ≡ 1 2 √ a 13 a 24 , where a ij = a i − a j , and where a 1 , a 2 , a 3 and a 4 are roots of the elliptic curve:

Analytic continuation
In this section we describe how to perform the analytic continuation to the physical region of linearly reducible elliptic Feynman integrals in a form that is suitable for fast and reliable numerical evaluations. Feynman integrals are analytically continued to the physical region by using the Feynman prescription, which is implemented by shifting the external invariants by a vanishing positive imaginary part iδ. Our starting point will be the one-fold integral representation of eq. (3.1) obtained from direct integration. Moreover we will assume that the integrand, at every ǫ order, is a pure polylogarithmic function of fixed transcendental weight multiplied by an overall algebraic function, of the form: , 1 x , where w 0 is the transcendental weight at leading ǫ order. In the equation above we made explicit the dependence on the Feynman prescription iδ which removes the branch cut ambiguities of the integrand in the physical region. y 2 is a quartic polynomial of x defining the relevant elliptic curve. The integrand of eq. (5.1) is symmetric under y → −y. This is due to the fact that the square root y appears when performing the integration with respect to the second-last Feynman parameter, by factorizing a certain second degree polynomial (this can be seen for example by analyzing the associated compatibility graph [13]), and the two roots of the polynomial are indeed symmetric upon flipping the sign of y. As we will show in the next sections the first master integral of the unequal masses sunrise topology (section 6.1) and the triangle with bubble integral (6.2) belong to this class. The analytic continuation of the second master integral of the sunrise topology can be done using the same techniques as described below. However, further analysis is required to obtain numerically stable representations due to the presence of simple poles in the leftmost integration kernels and we leave it for future work. Our task will be to identify a set of regions in the x, {s} space and remove, in each region, the dependence on the Feynman regulator δ by explicitly performing the δ → 0 limit: In order to perform the limits above we first need to compute lim δ→0 y(x, {s}, iδ). The square root y(x, {s}, iδ) for fixed {s} can be seen as a multivalued complex function of x and (vanishing) δ, taking two values differing by an overall sign. When defining the analytic continuation one usually defines a single-valued continuous branch of the square root for every x and δ, minus branch cuts for δ = 0 and a set of intervals x ⊂ R (see e.g. the discussion of [40]). There are multiple branches satisfying the constraints above, and one will have in general: where the actual signs depend on the definition of the branch under consideration. However, as discussed above, the integrand of eq. (5.1) is symmetric under y → −y and, in this case, the sign of the square root is immaterial and we set: The prescription above defines φ j (x, {s}) in every region R j and implies that in each region y 2 (x, {s}, 0) has definite sign.
In the examples discussed in the next sections we will consider polylogarithmic expressions up to weight three and, having a fast numerical evaluation in mind, we look for a representation of f j (x, {s}) can be found proceeding in the following algorithmic steps (see also [61,62]).
1. Function arguments are defined as monomials of the letters appearing in the symbol alphabet of f . In general also spurious letters might be needed when defining functions arguments (see for example [62]), and we have found this to be necessary for the upcoming examples. For the classical polylogarithms Li n (a(x, {s})), one requires that 1 − a(x, {s}) factorizes over the alphabet.
If the alphabet contains algebraic functions the factorization can be checked as follows. We take the logarithm of the function argument under consideration, and we equate it to an ansatz in the form of a linear combination of the logarithms of the alphabet letters. In this way we obtain a system of linear equations for the free coefficients of the ansatz. We numerically sample the equations for many values of the kinematic variables. If the equations admit a solution the argument factorizes as desired over the alphabet and the solution defines the factorized form.

JHEP01(2019)169
2. For each weight, one considers a set of linearly independent functions from the set of functions defined in the previous step. We have the freedom to choose the set of linearly independent functions defining the functional basis at weight i. We require that our basis elements satisfy: One then defines the most general ansatz for a Q-linear combination of these functions and products thereof, of weight i. The coefficients of the ansatz are then fixed imposing that the symbol of the ansatz equals the symbol of f (i) (x, {s}, 0).
3. We determine the terms in the kernel of the symbol at weight i by writing the most general ansatz in terms of the lower weight functions. We fix the free coefficients of the ansatz by specializing it to several points in the region under consideration. We then equate the resulting expressions to f (i) (x, {s}, iδ) and obtain a system of linear equations for the free coefficients that can be solved numerically with arbitrary precision. This gives a numerical value for the free coefficients of the ansatz that can be subsequently fitted against a basis of transcendental constants of appropriate weight.
The algorithm above does not rely on the rationality of the alphabet letters and generalizes the algorithm of [62] to algebraic cases.

Identifying admissible regions
In the previous section we have seen how to perform the δ → 0 limit at the integrand level and express the result in terms of analytic functions. The limit can be safely taken if we are able to identify a set of regions in the x, {s} space that contain no branch points for the integrand φ(x, {s}, 0)f (i+w 0 ) (x, {s}, 0). Let us show with an elementary example how a suitable set of regions can be identified. We consider the following elementary function: In order to be able to perform the limit we decompose the phase space in regions where the square root and the logarithm have no branch points: We can then explicitly perform the limit in the form: As we have seen in the previous section we will be interested in functions whose algebraic dependence comes only from the elliptic curve, therefore two regions are identified by requiring that the elliptic curve y 2 has definite sign:

JHEP01(2019)169
We then further partition these regions by requiring that the purely polylogarithmic expression, f (i) (x, {s}, 0), does not have branch points in the resulting subregions. The subregions are conveniently identified by studying the symbol alphabet letters. Specifically, the alphabet letters we will encounter for f (i) (x, {s}, 0) have the following general form: In general the partition above will overcount the number of regions that are actually needed. This is due to the fact that some of the zeros of the letters do not correspond to actual branch points of the polylogarithmic expression under consideration. While in principle one could perform a more refined analysis at this stage, for example by systematically studying the symbol map and the monodromy of the left-most symbol letters [63], in practice such overpartition is convenient when it comes to finding a basis of functions satisfying the constraints of eq. (5.7) in a given region. Indeed, it is true in general that the 'larger' the region the fewer are the functions satisfying the desired properties in that region, and in complicated cases one might find that the set of admissible functions does not span the functional space under consideration.

JHEP01(2019)169 6 Examples
In this section we provide a few examples that showcase the techniques discussed in the previous sections. We consider the unequal mass sunrise topology, and a triangle with bubble topology from both the direct integration and the differential equation point of view. By default we give our results in the Euclidean region. We provide the analytic continuation of the first master integral of the unequal mass sunrise, and of the triangle with bubble integral at order ǫ 0 and ǫ 1 in sections 6.1.3 and 6.2.3.

The off-shell sunrise diagram with unequal masses
We begin our discussion with the direct integration of the first master integral of the massive off-shell elliptic sunrise diagram with three different internal masses: where the subscripts of S 1,1,1 denote the powers of the internal propagators, and where we let s = p 2 . The same notation will be used in the remainder of the paper.

Direct integration
The Feynman parametrization of the sunrise has Symanzik polynomials: and is given by: We consider the case d = 2 − 2ǫ, where S 1,1,1 is finite. (Also note that the dimensionally regulated divergent 4-dimensional integral can be obtained from direct integration by analytic regularization [60], or from a dimensional recurrence relation.) In the Euclidean region p 2 < 0 and the internal masses are positive real valued. Expanding the integrand around d = 2 − 2ǫ gives: First we apply Cheng-Wu to put α 1 to 1, and integrate with respect to α 2 . The U polynomial is linear in the integration variable, whereas the F polynomial is not. Therefore we factor: where the roots are: (6.6) and we have a fourth degree polynomial in the last integration parameter: which defines an elliptic curve. Explicitly integrating with respect to α 2 one finds: At order ǫ 1 we obtain: where we introduced: In deriving this result we combined some logarithmic terms encountered at an intermediate stage.
Note that the inner polylogarithms are of uniform weight 2. Using HyperInt one may also obtain higher orders in ǫ with little difficulty. One can then verify by explicit computation that at order k in ǫ the polylogarithmic part is of weight k + 1.

Differential equations for the inner polylogarithmic part
Next we employ the ideas of section 4 to setup a system of differential equations for the inner polylogarithmic part. First we define the propagators explicitly: We may then write using the Feynman trick that: (6.12) LettingD 1 ≡ D 1 + xD 3 , we define a "generalized" topology where x is interpreted as an external scale: (6.13)

JHEP01(2019)169
which satisfies that: To perform an IBP reduction of the integrals in the inner polylogarithmic part we extend the topology with additional propagators to {D 1 , D 2 , N 1 , N 2 , N 3 }, where: and we obtain the IBP reduction using the C++ version of FIRE5. One may verify that in d = 2 − 2ǫ the following master integrals form a canonical basis: is the polynomial of eq. (6.7) with α 3 replaced by x. We have introduced a constant normalization for the inner polylogarithmic part by defining: (6.17) We note that we included the prefactor (m 2 3 ) 2ǫ in the canonical basis integrals to make them dimensionless. We divided out the term m 4 3 in the elliptic curve to obtain the form: where the a i variables denote the roots of the elliptic curve. In principle these are defined up to permutations, and for our purposes the ordering will not play an important role. To fix some convention for the roots, we let: , the canonical form differential equation is given by:

JHEP01(2019)169
where the differential equation matrix is: and where the letters l i are given by: We may obtain the symbol of the master integrals as long as we have their leading coefficients in the ǫ-expansion. One may find by power counting that B 3 vanishes at finite order. The leading coefficients of B 1 and B 2 are exactly 1. Therefore the leading coefficient vector is given by B (0) = (1, 1, 0). The symbol at all orders in ǫ can thus be written as: One may explicitly verify that the resulting symbol matches the symbol obtained from applying the maximal iteration of the coproduct to the solutions from HyperInt.
Solution in terms of multiple elliptic polylogarithms. We consider another method next, and solve the differential equation in terms of elliptic polylogarithms (E 4 -functions.) To do so we first map the integration parameter x to the domain [0, 1]. Note that for the initial application of Feynman's trick we could have alternatively used the form: We defineD 1 ≡ (1 − x ′ )D 1 + x ′ D 3 , and consider: and we may rewrite the canonical basis in terms of x ′ andŜ IPP as:

JHEP01(2019)169
where we now have the elliptic curve: The explicit expressions for the roots a ′ i are long and not particularly insightful expressions, and the reader may obtain them from the relation: We point out that the upcoming expressions in terms of E 4 -functions will be provided in the Euclidean region. This means we will use the following kinds of simplifications: We will solve the differential equation with respect to x ′ , which is given by: The partial derivative of A with respect to x ′ works out to: All of the entries may be expressed in terms of integration kernels of the E 4 -functions defined in [45], of which we wrote down the relevant ones down in eq. (4.18). We may write down the formal solution of the differential equation as a path-ordered exponential: and we find a particularly simple expression for the first master integral of the sunrise:

JHEP01(2019)169
which exposes the last integration kernel at all orders in ǫ. In order to obtain a representation in terms of E 4 -functions, we would like to pick the boundary condition x ′ 0 = 0, but we note that B(x ′ 0 , s, m 1 , m 2 , m 3 ) is singular in this limit. Nonetheless, the limit as x ′ 0 → 0 of the right hand side of eq. (6.33) should be finite, since the left hand side of the equation is finite.
Since the iterated integrals arising from the path-ordered exponential are multiple elliptic polylogarithms, we know that we may regulate the base-point divergence, which is of a logarithmic kind, using the tangential basepoint prescription. To get a consistent finite result we should apply the exact same regularization to the boundary term B(x ′ 0 , s, m 1 , m 2 , m 3 ), which will amount to taking the limit as x ′ 0 → 0 from the positive real axis, and throwing away divergences of the form log(x ′ 0 ) k , where k is a positive integer. Let's explicitly compute reglim x ′ →0 B(x ′ , s, m 1 , m 2 , m 3 ). It is relatively easy to compute the corresponding expression for B 1 , as the Feynman parametrization ofŜ IPP 2,0 has no non-trivial integrations. Furthermore, note that B 1 does not depend on s. For the integrals B 2 and B 3 there is a non-trivial integration to be performed. To compute the regularized limits of B 2 and B 3 we first exploit a symmetry based argument to simplify this integration.
If we put x ′ = 0 in the momentum space representation we find that the topology becomes that of a squared tadpole, and the resulting integral is independent of s. However, we need to first compute the integral for nonzero x ′ , and then compute the regularized limit as x ′ → 0 in order to get the correct boundary term. One may wonder if the dependence on s also disappears in the regularized limit, so that: We may write a closed form expression for B(x ′ , 0, m 1 , m 2 , m 3 ) by integrating up the Feynman parametrization. This leads to the following expressions: where we labeled the following terms: (6.37) Using these results we have tested numerically whether eq. (6.35) is justified for B 2 and B 3 . We took a random point in the external scales, with positive masses and a negative value for s. Furthermore, we took a number of increasingly small samples for x ′ , up to 10 −20 . We then computed B k (x ′ , s, m 1 , m 2 , m 3 ) for k = 2, 3 order by order in ǫ up to O(ǫ 3 ) for each value of x ′ by numerical integration. We note that in doing so it is important to first perform analytic regularization [60] of the Feynman parametrization so that the numerical integral converges order by order in ǫ. We then computed numerical values for B k (x ′ , 0, m 1 , m 2 , m 3 ) from the expressions in eq. (6.36). One finds that the difference B k (x ′ , s, m 1 , m 2 , m 3 ) − B k (x ′ , 0, m 1 , m 2 , m 3 ) becomes increasingly small for increasingly small x ′ (while the individual terms actually blow up as x ′ decreases.) By repeating this analysis for a few more points in the external scales we believe that eq. (6.35) is correct. Next we take the regularized limit as x ′ → 0. There are terms of the form (x ′ ) ǫ , which we first expand in ǫ, so that (x ′ ) ǫ = 1 + ǫ log(x ′ ) + 1 2 ǫ 2 log 2 (x ′ ) + O(ǫ 3 ), and then we throw away the logarithmic divergences. In other words we let the terms (x ′ ) ǫ → 1. The final expressions for the regularized limits are very simple pure functions of uniform transcendental weight: From eqs. (6.33), (6.34) and (6.38) we have all the elements to express B and in particular S 1,1,1 in terms of multiple elliptic polylogarithms. Note that in accordance with [45] we shuffle-regulate the E 4 -functions that end with a kernel of the type ψ 1 (0, x) = 1/x. For example we may write: where we explicitly worked out the shuffle product in one of the terms and used that E 4 ( 1 0 ; 1) = log(1) = 0. In terms of E 4 -functions the solution of the unequal mass sunrise in the Euclidean region is given up to order O(ǫ) 2 by: Lastly we remark on the three other master integrals in the top sector of the sunrise. One has for example: We can furthermore expressŜ IPP 2,2 in terms of the canonical basis integrals by IBP reduction,

JHEP01(2019)169
which leads to the following relation: where we have the following coefficients: and where the other coefficients are given by cyclic permutations: , for i = 1, 2, 3 and j = 2, 3, 4, and where we let a ′ 5 refer to a ′ 1 . It is clear from eqs. (6.33), (6.38) and (6.42) that S 1,2,1 can be integrated in terms of E 4 -functions. The first integrations are all expressible in terms of the kernels in eq. (4.18), and while the last integration contains kernels of the type dx ′ /(y ′ (x ′ − a ′ i )), it may be written in terms of kernels of E 4 -functions by IBP relations [45].

Analytic continuation
In this section we perform the analytic continuation to the physical region s > 0, m 2 1 > m 2 2 > m 2 3 > 0 of the first sunrise master integral S 1,1,1 (s, m 2 1 , m 2 2 , m 2 3 ) using the methods introduced in section 5. The analytic continuation of the ǫ 0 order is elementary as it requires only elementary identities among logarithms and we do not discuss it here, while we provide the result in appendix A. We discuss the analytic continuation of the order ǫ 1 coefficient, eq. (6.9), which for the reader's convenience we write here in the following form: where the iδ is introduced by applying the Feynman prescription s → s + iδ. The symbol alphabet letters of f (2) (x, s, m 2 1 , m 2 2 , m 2 3 , 0) can be expressed in terms of the following linearly independent letters:

JHEP01(2019)169
The alphabet of f (2) (x, s, m 2 1 , m 2 2 , m 2 3 , 0) contains only 8 linearly independent letters, however, as discussed in section 5, spurious letters are in general needed when defining function arguments, and in the present case we needed the extended alphabet above to be able to find a representation for lim δ→0 f (2) (x, s, m 2 1 , m 2 2 , m 2 3 , iδ) in the different relevant regions in terms of classical polylogarithms and logarithms (see e.g. [64]).
As explained in section 5 two regions are identified by requiring that the algebraic function, the elliptic curve y(x) in our case, does not have branch points: We then notice that in the region A S neither of the symbol letters vanish, so no further partitioning of A S is required. On the other hand region B S is partitioned as follows: Note that, as prescribed in section 5, each subregion is defined by requiring that all the letters have define sign, however for some subregions the set of constraints have no intersection, and only the subregions above need to be considered in this case. For later convenience we rename the regions as:

JHEP01(2019)169
The expressions for f (2) S,i (x, s, m 2 1 , m 2 2 , m 2 3 ), i ∈ {1, 3, 4, 5} are provided in appendix A. Note that in regions B S,j all the expressions have explicit imaginary parts and all the logarithms and dilogarithms are real valued, while this is not the case for region A S where the functions are complex valued in general due to the presence of i −y 2 in the arguments.
As already pointed out, the prescription of section 5.1 usually leads to an overcounting of the regions. This redundancy can sometimes be avoided by using the same set of functions (logarithms and dilogarithms in this case) for multiple subregions, and verifying that the resulting expressions are the same in different subregions. However these functions must satisfy the constraints of eq. (5.7), and in complicated cases, as the one under consideration, it is difficult to find a set of functions that satisfy these constraints on multiple subregions. Nevertheless this was possible for the triangle with bubble integral discussed in the next section.
Let us stress that the analytic continuation eq. (6.49) is suitable for fast and precise numerical evaluations, for example we have: In order to validate our results we performed extensive numerical checks against the computer program FIESTA. The results are summarized in figure 1.

Triangle with bubble
We consider the below triangle diagram, with a massive bubble insertion, relevant for the two-loop QCD corrections to heavy quark pair production: This diagram has the massive sunrise as a subtopology (as seen from contracting the massless internal propagator.) Hence we expect that the diagram cannot be expressed using multiple polylogarithms. Indeed an explicit calculation confirms this. In order to make the diagram finite in 4 dimensions we have put a dot on one of the massive propagators of the bubble. We note that there is only 1 master integral in the top sector of the topology to which this diagram belongs.

Direct integration
We take the following convention for the propagators: Then we have in particular: The Feynman parametrization has the Symanzik polynomials: 55) and is given by:

JHEP01(2019)169
We'll work in the Euclidean region where s < 0 and m 2 > 0. Expanding the integrand of T 1,2,1,1 (s, m 2 ) around d = 4 − 2ǫ gives: At order ǫ 0 one obtains: where we found it convenient to apply Cheng-Wu to set α 2 = 1. We can integrate with respect to the massless propagator to obtain: The polynomial m 2 (α 1 + α 3 + 1)(α 1 α 3 + α 1 + α 3 ) − α 1 α 3 s does not factor linearly in either of the remaining integration parameters without introducing a square root containing the other integration variable. Its roots are special cases of those encountered for the sunrise, namely R ± (s, m 2 1 , m 2 2 , m 2 3 ) corresponds to eq. (6.6) with α 3 replaced by α 1 . Performing the integration on α 3 yields us with: where we used the shorthand notation R ± = R Lastly, at order ǫ 1 we can integrate along the same path and we obtain the following result:
The higher orders in ǫ may be obtained from the same integration sequence.

Differential equations for the inner polylogarithmic part
We combine 2 massive propagators and define: where we let t = −s/m 2 be a scale with zero mass dimension. That way: We then have in particular: and we note that T 1,2,1,1 = T 1,1,2,1 by the symmetry of the diagram. Nonetheless, T IPP 3,1,1 = T IPP 2,2,1 are different polylogarithmic expressions, as they represent different integration sequences of the same integral. We adopt the notation: where P (x) S corresponds to eq. (6.7) with α 3 replaced by x. A canonical basis for the IPP is given by: where the coefficients are: The resulting canonical form differential equation is then given by: and where the letters are: l 1 = log (t) , l 2 = log (x) , l 3 = log (y) , l 4 = log (t + 1) , l 5 = log (x + 1) , l 6 = log (t + x + 1) , l 7 = log x 2 + tx + x − y + 1 x 2 + tx + x + y + 1 , We show next how to solve the resulting differential equation in terms of E 4 -functions. We perform the change of variables x = x ′ /(x ′ − 1), and let:

JHEP01(2019)169
where we picked the following convention for the roots: in the physical region t < −9. The differential equation matrix is given by: where the non-zero entries are: We let:T IPP a 1 +a 2 ,a 3 ,a 4 ≡ .

JHEP01(2019)169
Furthermore, we may write the full solution of B as a path-ordered exponential: and combining this with eq. (6.78) yields: We are interested in finding the boundary term: so that we may express eq. (6.80) in terms of E 4 -functions. One may verify that the top sector integrals B 1 and B 2 contain the termsT IPP 3,1,1 andT IPP 2,2,1 with prefactors that are proportional to an overall factor x ′ . Furthermore, we have the relations: and since T 1,2,1,1 = T 1,1,2,1 is a finite integral, the integrands in eq. (6.82) should have integrable singularities at the point 0. Therefore we find that: which we also verified numerically. One may compute B 3 for arbitrary x ′ by integrating the Feynman parametrization and one may observe that the limit x ′ → 0 vanishes as well. The canonical basis integrals B 4 , B 5 and B 6 belong to the sunrise subsector and their regularized limit may be obtained in the same manner as was done for the unequal mass sunrise topology. The results are: We now have almost all the ingredients to write T 1,1,2,1 in terms of E 4 -functions but there is a complication. Looking at eq. (6.80), one notices the appearance of the kernel ψ 1 (1, x ′ ) in the last entry. This kernel yields E 4 -functions that are individually divergent. First we consider the "naive" solution at finite order, which still contains divergent pieces:

JHEP01(2019)169
We would like to deal with the divergent terms by using the shuffle product to remove every occurrence of the kernel ψ 1 (1, x ′ ) in the first entry. A complication is that the kernel ψ −1 (1, x ′ ) may then appear in the first entry, which also diverges at 1. We deal with problem in a similar manner to [45], where such issues arise in the analysis of the second master integral of the equal mass sunrise. First, we define a new kernel: which is a regulated version of ψ −1 (1, x ′ ). We then express our E 4 -functions in terms of this new kernel. After doing so one may extract out the divergent pieces from each E 4 -function by shuffle regularization. The only remaining divergent terms are: and their prefactors should vanish as we know T 1,1,2,1 is finite. One finds for example the contribution: dmath and it may be numerically verified up to high precision that the combination of E 4 -functions multiplying E 4 ( 11 11 , 1) evaluates to zero. We decide to restore the kernel ψ −1 (1, x) in all entries but the first, and we obtain the following representation in terms of E 4 -functions that are individually finite: The higher orders in ǫ may be obtained in the same manner.

Analytic continuation
In this section we apply the methods of section 5 to perform the analytic continuation of the triangle with bubble integral T 1,2,1,1 (s, m 2 ) to the physical region s > 0, m 2 > 0. We will explicitly discuss the analysis for the order ǫ 0 contribution T 1211 . We also performed the analytic continuation of the order ǫ 1 , however the final result involves rather JHEP01(2019)169 complicated expressions and we don't present its derivation here. Our starting point is the representation of eq. (6.60) that we rewrite here for the reader's convenience in the following form: T (x, s, m 2 , iδ)dx , (6.90) where we renamed the integration variable to x, and where f (2) (x, s, m 2 , iδ) is obtained from eq. (6.60) by applying the Feynman prescription s → s + iδ. Referring to section 5.1, the symbol alphabet letters of f (2) (x, s, m 2 , 0) can be expressed in terms of the following linearly independent letters (see section 6.1.3 for further discussion): Following the prescription of section 5.1, we identify the following subregions of region A T : y 2 (x, s, m 2 ) < 0: A T,1 : ρ 4 > 0 , A T,2 : ρ 4 < 0 , (6.92) while for B T : y 2 (x, s, m 2 ) > 0, we have the following subregions: By applying the algorithm of section 5 we are able to perform the δ → 0 limit of f (2) (x, s, m 2 , iδ) in terms of classical polylogarithms and logarithms in each of these regions. We obtain: 1 (x, s, m 2 ) is, by using the prescription of eq. (5.6), related to the one of f 3 (x, s, m 2 ) by: while: + log (ρ 2 ) log (ρ 5 ) + log(2) log (ρ 5 ) + log 2 (−σ 2 ) − 2 log(2) log (−σ 2 ) + log 2 (2) , (6.97)
T (x, s, m 2 , iδ) in each of these "enlarged" regions leads to the same expression for all the corresponding subregions, therefore only the 3 regions R T,i are needed.
Since x, s, y are real valued in R T,2 and R T,3 , and since by construction the logarithms and dilogarithms of eqs. (6.97), (6.98) satisfy eq. (5.7), the terms f (2) T,2 (x, s, m 2 ) and f (2) T,3 (x, s, m 2 ) have explicit imaginary parts and all the logarithms and dilogarithms are real valued. This is not the case for f (2) T,1 (x, s, m 2 ) where the dependence on i −y 2 implies that individual functions are complex valued in general.
Let us mention again that the analytic continuation eq. (6.94) is suitable for fast and precise numerical evaluations, for example we have:

A non-planar triangle from Higgs+Jet
In this section we show that linear reducibility does not directly imply that a representation in terms of eMPL exists, and further exploration of the methods discussed in this paper is required. Nevertheless we provide evidence that a simple all orders structure, analogous to the one discussed in the previous section, holds. We consider a non-planar triangle in d = 4 − 2ǫ with two off-shell legs, relevant for the two-loop QCD corrections to Higgs plus jet production: whose homogenous solutions were found in [37]. This integral has an elliptic maximal cut, but no elliptic subtopologies are present. Note the (1+2ǫ)(p 2 2 −s) prefactor, which is needed to obtain a uniformly transcendental inner polylogarithmic expression. We consider the Euclidean region s < 0, p 2 2 < 0 and m 2 > 0. The parametric representation in d = 4 − 2ǫ dimensions reads:

Conclusions and outlook
In this paper we have investigated a class of elliptic Feynman integrals, dubbed linearly reducible elliptic Feynman integrals. By direct integration of the Feynman parametrization one may express such integrals order by order in the dimensional regulator as a 1dimensional integral over a polylogarithmic integrand, which we call the inner polylogarithmic part (IPP). The resulting 1-dimensional integral representation can be analytically

JHEP01(2019)169
continued to the physical region in a form suitable for fast and precise numerical evaluations. When the IPP depends on one elliptic curve and no other algebraic functions, linearly reducible elliptic Feynman integrals can also be expressed in terms of multiple elliptic polylogarithms.
We have also shown that the IPP can be mapped to a (generalized) polylogarithmic Feynman integral that can be subsequently solved using the differential equations method. In particular we studied the IPP of the unequal mass sunrise topology, and a triangle with bubble topology, and provided a canonical basis of master integrals where the system of differential equations is in d log ǫ-factorized form. For this basis, the differential equations with respect to the last integration parameter were found to be in an ǫ-form where the integration kernels coincide with the integration kernels of the class of eMPLs of [44]. This allows one to systematically solve the IPP in terms of eMPLs, directly from the system of differential equations. Once such a representation is achieved the remaining last integration can be performed in terms of eMPLs as well.
We expect that the methods discussed in this paper may provide new insights for problems where iterated integrals over multiple (and more complicated) algebraic functions need to be considered. Furthermore we aim to apply our methods to more complicated Feynman integrals in the future.