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 $\epsilon$-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.


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. Furthermore, Feynman diagrams have exposed a rich mathematical structure. Most progress has come from identifying a class of Feynman integrals expressible in terms of special functions known as multiple polylogarithms (MPL) [1,2], to which all 1-loop scalar diagrams, all 2loop scalar massless on-shell 4-point diagrams, and many more belong [3]. 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 [4] and integration-by-parts (IBP) identities [5][6][7][8] 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.
From the traditional viewpoint Feynman integrals are defined as momentum space integrals: one integrates over a 4-dimensional space for each internal momentum. Furthermore, in dimensional regularization the space is taken to be d-dimensional. Unfortunately from this viewpoint Feynman integrals are difficult to calculate both analytically and numerically. To get a handle on Feynman integrals it is therefore necessary to provide simpler representations of them. Arguably the two most successful but different starting points are the differential equations method [9][10][11][12][13] and the direct integration of the scalar parametrization of Feynman integrals [3,14].
The differential equation method has received significant attention in the last decade. In this method collections of Feynman integrals are considered which are taken to be closed under differentiation with respect to external scales, like the Mandelstam variables and the masses. The closure of such sets under differentiation is shown by employing IBP identities. Solving the resulting differential equations in the case of multiple polylogarithmic integrals has been systematized by the introduction of the so-called canonical basis [15]. In the polylogarithmic case the main challenge is reducing the generally very large system of IBP equations that is necessary to obtain a closed form for the differential equation. Specialized software packages are available to perform these reductions [16][17][18], but all run into computational limits for sufficiently complicated topologies. The problem of finding the canonical basis has been understood quite well in the meantime and has been automated in a number of software packages [19][20][21]. On the other hand the direct integration of the scalar parametrization of Feynman integrals has seen major advancements by the works of (among others) Brown and Panzer, see for example [3,14]. In particular the Maple package HyperInt [22] automates most of the direct integration methods, and we use it to perform the integrations in this paper.
At 1-loop a (multiple) polylogarithmic basis is formed by the bubble, triangle, box and pentagon scalar diagrams [3,23]. At 2-loop and higher there is no canonical way to define a basis, and diagrams are not in general expressible in terms of multiple polylogarithms. Inclusion of a diagram in the class of multiple polylogarithmic diagrams can be shown, e.g., by proving that the compatibility graph associated to the parametric representation is linearly reducible [3,14]. The next challenge is to get a grasp on the analytical structure of non-MPL integrals. The first step is to understand the class of elliptic Feynman diagrams . Their name follows due to the fact that their maximal cut is expressible in terms of combinations of elliptic integrals. These integrals furthermore satisfy a second or third order differential equation. More complicated classes of Feynman diagrams will satisfy even higher order differential equations and are a challenge left for the future.
Considerable progress has been made in the understanding of elliptic Feynman diagrams. The general strategy is to define an elliptic generalization of multiple polylogarithms. Quite different definitions of these functions have been made in the literature [30,[45][46][47]. Elliptic multiple polylogarithms are defined by requiring that the integration kernels are independent with respect to integration by parts identities and that they have only simple poles. Moreover these functions are sufficient to express all iterated integrals over an elliptic curve. In the context of the differential equations method a recent insight has proven useful. Because 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). After solving for the subtopologies one can obtain an expression for the full integral as a simple application of the variation of parameters method. The computation of maximally cut Feynman integrals has been aided by recent explorations of the Baikov parametrization [48][49][50].
Nonetheless a complete understanding of elliptic Feynman integrals is still missing. In this paper we investigate dimensionally regulated elliptic Feynman integrals by direct integration of the Feynman parametric representation. In so doing we identify a class of so-called linearly reducible elliptic Feynman integrals, that can be algorithmically solved up to arbitrary order of the dimensional regulator in terms of a suitable elliptic generalization of multiple polylogarithms. More specifically, the direct integration approach requires that the polynomials defining the Feynman parametric representation factor linearly with respect to the integration variable (Feynman parameter). However, in the elliptic case, the factorization introduces irreducible square roots 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. We then argue that linearly reducible elliptic Feynman integrals can be expressed, to all orders, in terms of the following elliptic generalization of multiple polylogarithms, . . , a n (α); 1) , where n is an integer, a 0 is a constant with respect to α, δ ∈ {0, 1}, G (a 1 (α), . . . , a n (α); 1) is a multiple polylogarithm of argument 1, a 1 (α), . . . , a n (α) are algebraic functions in α, and P (α) is a polynomial in α. The definition of a (Goncharov) multiple polylogarithm G is given in eq. (2.6). While on general grounds the functions (1.1) are expected, we found that for all the examples discussed in the main text only the special cases where n ∈ {0, 1} and P (α) is a polynomial of degree 3 or 4 with no repeated roots, encoding the information of the elliptic curves associated to the integral, are needed. Moreover, for the upcoming examples it was always possible to perform an integral basis choice such that the results at a given ǫ order are of the general form, where UT(α p ) is a pure polylogarithmic function of uniform weight and P el is a polynomial of degree 3 or 4 with no repeated roots. While we do not have a proof that this property holds to all orders in ǫ, we have verified it up to high orders. 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 define linearly reducible elliptic Feynman integrals, and discuss the class of functions to express such integrals. In section 4 we consider multiple examples that all fall into this class of functions. In appendix A we solve a next-to-linearly reducible integral, and we show that a two dimensional generalization of (1.1) is required in this case. In appendix B we provide explicit full results for the most complicated examples discussed in the main text. We end our discussion in section 5, where we reflect on the direction integration method and our results and compare them to other available methods. Lastly we provide an outlook on future work.

Parametric representation of Feynman integrals
The defining momentum space integral representation of scalar Feynman diagrams admits the form: where {s} = {p} ∪ {m} 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 on the propagators are integers.
A scalar parametrization can be found as follows. First one uses the Schwinger trick to rewrite every propagator as: 1 1 The momentum integrals of eq. (2.1) become Gaussian integrals which can be performed to give: 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 2 F 2 (G) of the Feynman diagram G: where s (T1,T2) 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 [51,52]. 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 choose the normalization constant to remove the prefactor in the Feynman representation.

Direct integration, linear reducibility and all orders statements
First we remind the reader of the definition of (Goncharov) multiple polylogarithms. They are the following recursively defined functions: 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. We have the definition: to deal with the divergent case where all n parameters a i are equal to zero. The recursive definition in eq. (2.6) comes to a halt for n = 0, where we let in agreement with eq. (2.7): A pedagogical review of these functions and their functional identities is given in [53]. The direct integration of the parametric representation discussed in the previous section can be performed, e.g., by using the computer program HyperInt [22]. In general, the parametric representation requires, order-by-order in ǫ, the computation of p integrals of the form, 1 Convergence of the scalar integral is guaranteed by adding a small imaginary part to the propagators according to the Feynman prescription. We'll ignore this subtlety in the remaining text.
2 A 2-forest is a disjoint union of 2 trees.
where k ∈ {1, . . . , p} and f k (α k+1 , . . . , α p ) is a combination of multiple polylogarithms, which depends on a set of irreducible polynomials that we denote by P (k) (α k+1 , . . . , α p ). Schematically, the integration of (2.9) involves the following three steps (see refs. [22] and [54] for further details), and requires that P (k−1) (α k , . . . , α p ) depends in a rational way on the integration variable α k , 1. express f k−1 (α k , . . . , α p ) as a combination of multiple polylogarithms of argument α k , The reason why P (k−1) (α k , . . . , α p ) needs to be rational in α k is that steps 1 and 2 rely on its linear factorization with respect to α k , where we denoted by P (k−1) i the i-th polynomial of the vector P (k−1) , and by m i the degree of the polynomial P In some cases the g k−1 i,j (α k+1 , . . . , α p ) are algebraic functions (e.g. they contain square roots). This can be avoided sometimes by applying the Cheng-Wu theorem (see next section), choosing a different integration sequence, or performing a variable change that removes the algebraic functions, otherwise the algorithm stops and the next integration of f k (α k+1 , . . . , α p ) with respect to α k+1 in terms of multiple polylogarithms is not possible.
Let us now discuss more precisely the parametric integration in the case of elliptic Feynman integrals. By definition, an elliptic Feynman integral depends on one or more integrals over a set of elliptic curves, i.e. irreducible square roots of polynomials of degree 3 or 4 with no repeated roots. It is therefore clear that in the elliptic case there exists an α k such that g k−1 i,j (α k+1 , . . . , α p ) depends on a set of irreducible square roots, and the remaining integrations with respect to α k+1 , . . . , α p are not possible. We define linearly reducible elliptic Feynman integral an integral such that, to all ǫ orders, and where the algebraic functions depend on at most square roots (i.e. no cubic or higher roots are allowed). Therefore the final expression of a linearly reducible elliptic Feynman integral reads, to a generic ǫ order, where f k (α p ) 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. Nonetheless this definition is motivated by the fact that we find empirically that many elliptic Feynman integrals are indeed linearly reducible. A remarkable feature of the direct parametric integration is that it is easy to make all orders statements by means of the compatibility graph method introduced in [14]. While we refer the reader to that reference for further details, the main consequence of the compatibility graph method is that the polynomials P (k) (α k+1 , . . . , α p ), k ∈ {1, . . . , p} are independent of the ǫ order we are considering, starting from the next-to-leading order of ǫ. This implies that if, e.g., a finite elliptic Feynman integral is linearly reducible at order ǫ 1 , it is linearly reducible to all orders. Moreover, the alphabet of the polylogarithmic expression f k (α p ), is the same to all orders. In section 4 we solve a number of finite linearly reducible elliptic Feynman integrals up to and including the order ǫ 1 , and it is understood that this proves linear reducibility up to all orders.

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 is a nonempty set of integers [55]. Projectivity refers in this context to the property that the integrand should remain invariant under the rescaling: If one starts with an integral over ∆ which does not remain invariant under eq. (2.14), the following projective transform can be performed to obtain such an integral [56]. 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: (2.16) Eq. (2.16) 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. 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: 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, labelled by a set S ′ , and where the integrand is projective. By the Cheng-Wu theorem we can turn this into:

The class of functions
In the previous section we argued that linearly reducible elliptic Feynman integrals take, to all orders ǫ, the schematic form: where the sum over i denotes a generic collection of terms grouped by factors A i which are algebraic functions in α, and MPL i (α) is a polylogarithm with algebraic arguments. As defined in the previous section the algebraic functions contain at most square roots, and no higher roots are allowed. Using the compatibility graph method, elementary identities between algebraic functions and partial fractioning, one can show that (3.1) can be expressed as a linear combination of the following elliptic generalization of multiple polylogarithms, a 1 (α), . . . , a n (α); 1) , where n is an integer, δ ∈ {0, 1}, a 1 (α), . . . , a n (α) are algebraic functions in α, and a 0 is a constant with respect to α. While the functions (3.2) are expected to be required to solve a generic linearly reducible elliptic Feynman integral, it is remarkable that for the examples discussed in the next section only the special cases with n ∈ {0, 1} and where P (α) is a polynomial of degree 3 or 4 were needed, by performing a proper choice of powers of propagators and of dimension. Moreover it will be interesting to investigate the possibility of systematically reducing the general case (3.2) to a minimal set of special cases by using partial fractioning and integration-by-parts identities as done in [57] for the generalized harmonic polylogarithms. We find it then natural to define: a 1 (α), . . . , a n (α); 1) , where a 1 (α), . . . , a n (α) are algebraic functions in α and P el (α), and a 0 is a constant with respect to α. The superscript α is included to remind ourselves of the integration parameter in the arguments. Special cases of this class of function have been found recently in [41,42] for the imaginary part of the sunrise integral and the 10 points elliptic double box. The previous definition also covers cases in which the elliptic curve or (α − a 0 ) are absent in the algebraic integration factor, and in these cases we use respectively the shorthands: G α ( P el (α); a 1 (α), . . . , a n (α)) ≡ ∞ 0 dα 1 P el (α) G (a 1 (α), . . . , a n (α); 1) , G α (a 0 ; a 1 (α), . . . , a n (α)) ≡ G (a 1 (α), . . . , a n (α); 1) .

(3.4)
In fact, all the linearly reducible examples discussed in section 4 were expressible in terms of the two special cases (3.4), while the general case (3.3) was needed for the next-to linearly reducible example discussed in appendix B, along with a suitable two dimensional generalization thereof. Our examples sometimes involve quite lengthy expressions, so we use the following compact notation for multiple polylogarithms evaluated at 1 and their elliptic generalization (3.3), G a1,...,an ≡ G(a 1 , . . . , a n ; 1), G α √ P el ,a0;a1(α),...,an(α) ≡G α ( P el , a 0 ; a 1 (α), . . . , a n (α)), (3.5) which makes full results suitable to be presented explicitly.
The numerical evaluation of functions in the class of eq. (3.3) can be performed by combining a numerical integrator with a library, e.g., GiNaC, that is capable of fast evaluation of the inner polylogarithms. We have used this method to numerically cross-check our results with FIESTA [58].

Examples
In this section we provide a number of examples where we apply the techniques discussed in the previous sections. All the integrals are solved up to and including at least the order ǫ 1 , which proves linear reducibility to all orders. Moreover, the computation of the higher orders is completely algorithmic and can be done by applying the Cheng-Wu theorem and the integration order used for the order ǫ 1 . Remarkably, linear reducibility was always achieved combining the Cheng-Wu theorem with a suitable integration order, and no variable change was necessary -in fact, finding a proper variable change is in general a very non-trivial task, while finding the right integration order and Cheng-Wu replacement can be automated. In fact we checked that many more examples than the ones we present in this paper do indeed belong to the class of linearly reducible elliptic Feynman integrals, e.g. the kite integral and the Higgs+Jet non planar elliptic topologies. All the results are provided for the Euclidean region where all the external invariants are negative real valued and the internal masses are positive real valued. To avoid clutter we set the gamma function and the π normalization factors of the Feynman parametrization to 1, by the choice of normalization in eq. (2.5).

The sunrise diagram
We start our discussion with the simplest elliptic Feynman integral, the massive off-shell sunrise with three different internal masses: where the subscripts of S 111 denote the powers of the internal propagators. The same notation will be used in the remainder of the paper. The Feynman parametrization of the sunrise has Symanzik polynomials: and is given by We consider the case d = 2 − 2ǫ, where the sunrise is finite. (Note that the dimensionally regulated divergent 4-dimensional integral can be obtained from direct integration by analytic regularization [59], or from dimensional recurrence relations.) Expanding the integrand around d = 2 − 2ǫ gives: We apply Cheng-Wu to put α 3 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: and we have a fourth degree polynomial: which is an elliptic curve in the last integration variable. Explicitly integrating with respect to α 2 one finds: At order ǫ 1 we obtain: where we introduced: (4.10) In deriving this result we combined the logarithmic terms encountered at an intermediate stage. Note that the inner polylogarithms are of uniform weight 2. Similarly one may verify by explicit computation that at order k in ǫ the polylogarithmic part is of weight k + 1. At order ǫ 2 one obtains the following result: where we introduced the term:

(4.18)
At this point it is straightforward to perform the last integration with respect to α 4 . In this way the final result is an integral with respect to α 5 between 0 and 1. In order to express the result in terms of the functions eq. (3.2) we perform the following variable change which maps the upper integration bound to ∞ as desired. The final expression for N (0) 111111 (s, p 2 2 , m 2 ) is quite lengthy and we provide it in appendix B. Its expression is, as expected, of the form: where the elliptic curve is The explicit result is very lengthy and can be obtained upon request from the authors.

Triangle with bubble
We consider the below triangle diagram, with a massive bubble insertion, relevant for the two-loop QCD corrections to the 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 put a dot on one of the massive propagators of the bubble. The Feynman parametrization has the Symanzik polynomials: Expanding the integrand 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 , m 2 , m 2 ), where R ± (p 2 , m 2 1 , m 2 2 , m 2 3 ) is defined in eq. (4.6). Performing the integration on α 3 yields us with: where we used the shorthand notation R ± = R ± (s, m 2 , m 2 , m 2 ), and have introduced the term: Lastly, at order ǫ 1 we can integrate along the same path and we obtain the following result: where we introduced the additional terms: . (4.32)

Conclusions and outlook
In this paper we have identified a class of elliptic Feynman integrals, dubbed linearly reducible elliptic Feynman integrals, that can be algorithmically computed to all orders of the dimensional regulator in terms of a suitable elliptic generalization of multiple polylogarithms, by directly integrating their Feynman parametric representation. In particular, the results for linearly reducible elliptic Feynman integrals can be expressed, to all orders, as a one-dimensional parametric integral over a polylogarithmic function (whose maximal transcendental weight increases one unit at each order) depending on the relevant elliptic curves via its arguments and coefficients. The all orders statements of the Feynman parametric representation are easily proven by using the compatibility graph method developed in [14].
For the integral topologies we considered, it was always possible to choose a representative integral and a normalization factor such that the result is the integral over a pure polylogarithmic function of uniform weight, multiplied by an overall algebraic factor. While this property was not proven to all orders, we have directly checked, for our examples, that it holds up to relatively high orders. It will be interesting to develop general methods to construct an integral basis satisfying such a generalized "canonical" form. In fact, for the integrals we considered an overall normalization factor was sufficient to achieve this property, while in general we expect to need linear combinations of integrals. We leave this investigation for future work.
The direct integration strategy studied in this paper has proven to be very efficient when compared, e.g., to the differential equations method, as it does not require the solution of (generally large) systems of IBP relations, and the computation of boundary conditions. Moreover, in the elliptic case, higher order (homogeneous) differential equations need to be solved, which is in general an unsolved problem. On the other hand in our direct integration approach the information about the relevant elliptic curves is naturally encoded by the roots of sets of polynomials encountered during the integration.