Analytic results for two-loop planar master integrals for Bhabha scattering

We analytically evaluate the master integrals for the second type of planar contributions to the massive two-loop Bhabha scattering in QED using differential equa- tions with canonical bases. We obtain results in terms of multiple polylogarithms for all the master integrals but one, for which we derive a compact result in terms of elliptic mul- tiple polylogarithms. As a byproduct, we also provide a compact analytic result in terms of elliptic multiple polylogarithms for an integral belonging to the first family of planar Bhabha integrals, whose computation in terms of polylogarithms was addressed previously in the literature.


Introduction
The Bhabha scattering process, i.e., the elastic scattering of an electron-positron pair, is one of the standard candles at lepton colliders, and it will play a crucial role at future circular or linear colliders. A precise theoretical knowledge of this process, including up to nextto-next-leading order effects in the QED coupling constant, is therefore highly desirable. So far the NNLO cross section is only known in the massless limit, supplemented by the leading logarithmic finite mass effects (see, e.g., ref. [1] for recent results). Complete NNLO results including the full dependence on the electron mass, however, are not yet available.
One of the main obstacles to obtain the complete NNLO results is the complexity of the two-loop integrals involved. The relevant two-loop integrals were evaluated in the smallmass limit in refs. [2][3][4][5][6]. Up to now, there are only partial results for two-loop integrals with massive fermions. Analytic results for diagrams with one-loop insertions and a closed massive fermion loop were obtained in refs. [7,8]. First analytic results for massive twoloop double-box Bhabha diagrams were obtained in refs. [9,10]. More attempts to evaluate two-loop Bhabha integrals can be found in refs. [11][12][13][14][15].
The master integrals relevant for the calculation of the two-loop corrections to Bhabha scattering can be classified into two planar families and one non-planar one. The systematic evaluation of the integrals in the first family was started in ref. [16], where the master integrals for the family associated with graph (a) of fig. 1 were evaluated in terms of multiple polylogarithms [17][18][19] (MPLs) in the framework of the method of differential equations [20][21][22] with the help of the strategy based on canonical bases [23]. In ref. [16], a solution in terms of MPLs had been provided for all integrals except one, whose analytic evaluation was hindered by the presence of the a non-rationalisable square-root in the symbol alphabet. More recently, it was shown that also this last integral can be expressed in terms of multiple polylogarithms by an integration technique based on an ansatz of MPLs with suitable arguments [24]. The goal of the present paper is to analytically evaluate the master integrals for the second planar family, which is associated with graph (b) of fig. 1. In a way that is reminiscent of the first planar family of integrals, we will show that also in this case we can obtain results in terms of MPLs for all master integrals but one (see fig. 2), due to the presence of the same non-rationalisable square root found in the evaluation of the master integrals for graph (a). While it can be shown by direct integration techniques that also for this master integral a representation in terms of MPLs exist, the representation we obtained is extremely cumbersome and of no practical use. Nevertheless, it turns out that a compact result for this integral can be derived in terms of the elliptic MPLs introduced in refs. [25][26][27]. As a byproduct of our analysis, we will also provide a very compact analytic result for the remaining master integral in the first family in fig. 1(a) in terms of the same class of functions. 1 The rest of the paper is organised as follows. In section 2 we explain the notation, present the system of canonical differential equations for the problem at hand, and introduce the alphabet needed to describe the planar family in fig. 1(b). The alphabet contains multiple square roots, and it is a priori not obvious that the differential equations can be integrated in terms of multiple polylogarithms. In section 3 we review general sufficient criteria and algorithms to solve a system of canonical differential equations in terms of MPLs or their elliptic generalisations. We then continue in section 4 to apply these ideas to our computation. We explain how our differential equations can straightforwardly be solved in terms of multiple polylogarithms for all elements of the basis of the master integrals but one, which is connected with the graph shown in fig. 2. In section 5 we show how a convenient and compact representation for this integral can be found in terms of elliptic generalisations of MPLs. We also provide a very compact, alternative, analytic solution for one of the master integrals in the first planar family considered in refs. [23,24]. Finally we draw our conclusions in section 6.

Canonical differential equations
The Feynman integrals for the family of fig. 1 (b) can be organised in an integral family with nine propagators, where the first seven propagators correspond to the edges of the graph and the last two are so-called irreducible numerators, (2.1) The indices a i , i = 1 . . . 7 can be positive or negative integers, but we restrict our computation to a 8 , a 9 ≤ 0. We work in dimensional regularisation in D = 4 − 2 dimensions in order to regulate both infrared and ultraviolet divergences. The electron mass is denoted by m, and the external momenta p i are on shell, p 2 i = m 2 . We introduce the usual Mandelstam variables Using the public codes FIRE [28] and KIRA [29,30], we solve the integration-by-parts (IBP) identities [31,32] and reveal 43 independent master integrals, which we collect into the vector g = (g 1 , . . . , g 43 ) T . This vector satisfies a system of linear differential equations of the form where v = s, t, m 2 , ∂ v = ∂ ∂v and the matrices A s , A t , A m 2 are rational functions of s, t, m 2 and . In the following, it will be useful to collect all the partial derivatives into a total differential, and to work with the equation In order to evaluate the master integrals, it is convenient to search for a so-called canonical basis [23], i.e., a basis of master integrals for which the matrix dA on the righthand side of eq. (2.4) is proportional to and its entries can all be expressed as linear combinations of total differentials of logarithms. While it has been suggested (at least conjecturally) that the study of the residues of the integrands can provide the full information to determine the elements of the canonical basis [33], this analysis can become computational very expensive when square-roots are involved. 2 Therefore, we follow here a mixed approach to find a canonical basis. In order to select suitable candidates, we start by analysing only the leading singularities associated to the maximal cuts of the various integrals, and we supplement this analysis by the method of ref. [35]. By choosing integrals with unit leading singularities at the level of the maximal cuts, one can often bring the initial differential equations into a so-called precanonical form, where the corresponding matrices depend linearly on . Once this is achieved, the prescriptions of ref. [35] can be successfully applied to arrive at a fully canonical form. In our case, the precanonical basis is composed of the 43 Feynman integrals present on the right-hand side of the above expression for the -basis. In this way we obtain the new system of differential equations df = dĀ f , (2.5) withĀ independent of . Our canonical basis f i is given in appendix A. The price to pay for casting the equations in canonical form is that the new matrix dĀ is not a matrix of rational one-forms, but it involves four square roots: More precisely, the matrixĀ takes the general form where the R i are algebraic functions, referred to as letters (and their collection is called 2 See ref. [34] for an automated implementation of this approach. the alphabet): As we will see, the appearance of these square roots makes the solution of the differential equation (2.5) in terms of MPLs highly non-trivial. It is easy to write down a formal solution to eq. (2.5) as where Pexp denotes the path-ordered exponential and f 0 ( ) encodes the initial condition and is related to the value of f at a specific point (s 0 , t 0 , m 2 0 ). The path γ connects the initial point (s 0 , t 0 , m 2 0 ) to the generic point (s, t, m 2 ). When expanding eq. (2.9) in , then at each order we can write f in terms of Chen iterated integrals [36], defined in the following way: consider a path γ and a collection of one-forms ω i . If t denotes a local coordinate on γ, we can pull each ω i back to γ and write γ * ω i = dt F i (t). We can then define the iterated integral of a sequence ω i 1 . . . ω in (usually referred to as a word ) by (2.10) The number n of integrations is called the length of the iterated integral. In general, this integral will not be homotopy-invariant, i.e., it will depend on the details of the path γ. There is a necessary and sufficient condition, called the integrability condition, for a combination of iterated integrals to be homotopy-invariant [36]. The details of this criterion are not important in the following. Here it suffices to say that it is always satisfied for the solutions in eq. (2.9). The corresponding Chen iterated integrals are then multi-valued functions of the end point (s, t, m 2 ) of the path γ, where the multi-valuedness only comes from choosing two non-homotopic paths from (s 0 , t 0 , m 2 0 ) to (s, t, m 2 ). The master integrals f (s, t, m 2 ; ) can be expressed at every order in terms of Chen iterated integrals over words from the alphabet in eq. (2.8). The coefficient of n in the path-ordered exponential in eq. (2.9) only involves iterated integrals of length n. Chen iterated integrals are a very general class of functions, and it can be useful to express them in terms of a class of functions that are well-studied in the literature. In the next section we discuss sufficient criteria for when Chen iterated integrals over d log-forms with algebraic arguments, as the ones above, can be expressed in terms of other classes of functions.

From d log-forms to multiple polylogarithms
Since the solution of differential equations in canonical form as Laurent series in dimensional regularisation leads naturally to linear combinations of Chen iterated integrals over words of d log-forms, it is natural to ask when it possible to evaluate such iterated integrals in terms of other classes of functions, and if it is possible to perform this rewriting algorithmically. The advantage of this rewriting lies in the fact that these classes of functions may be well studied in the literature, and there may be established methods or computer codes for their manipulation and/or their numerical evaluation.
Let us consider Chen iterated integrals of the form γ w, where γ is a path from an initial point x 0 = (x 0,1 , . . . , x 0,p ) to the point x 1 = (x 1,1 , . . . , x 1,p ), and w is an integrable combination of words of length n of the form We consider the initial point x 0 fixed, and we see the integral as a multi-valued function of the end-point x 1 . The arguments of the logarithms are assumed to be algebraic functions of the variables x = (x 1 , . . . , x p ). For some time, there was a folklore belief in the physics community that all such iterated integrals could be evaluated in terms of a rather simple class of iterated integrals, namely the so-called multiple polylogarithms (MPLs), defined by [17][18][19]37] G a 1 ,...,an (z) = G(a 1 , . . . , a n ; z) = z 0 dt t − a 1 G(a 2 , . . . , a n ; t) , where the recursion starts with G(; z) ≡ 1. In the special case where all the a i 's are zero, we define G 0,...,0 (z) = G(0, . . . , 0 The shared belief in physics was that, whenever the d log's in eq. (3.1) have algebraic arguments R i , then this integral can be written in terms of MPLs whose arguments a i and z are algebraic functions of x 0 and x 1 . This belief was shown to be false in ref. [38], where an example of an iterated integral of length two was constructed that cannot be evaluated in terms of MPLs with algebraic arguments. The result of ref. [38] shows that the question of whether an iterated integral of d log-forms with algebraic arguments can be expressed in terms of MPLs depends in general on the details of the integral, including the details of the integration contour.
In the remainder of this section we discuss some special cases of algebraic letters R i where the Chen iterated integrals can be evaluated in terms of other classes of special functions. The starting point is the observation that in the context of Feynman integrals, the R i are not generic algebraic functions, but often involve at most square roots of polynomials, i.e., they are of the form: where the r i,j are rational functions, and the q a are polynomials. Without loss of generality, we assume the r i,j to be polynomials. At this point we have to make a comment: The integrand of the iterated integration, and thus the functional form of the letters R i , is sensitive to a change of variables. In particular, one may ask if we can find a change of variables x = ψ(y), with ψ a rational function, such that q a (ψ(y)) = u a (y) 2 is a perfect square, for some values of a. This operation is known as rationalisation of square roots. Over the last few years, several algorithmic criteria have been developed to find such a parametrisation for specific Feynman integral computations, or to prove that rationalisation is instead not possible [39][40][41][42][43][44][45]. Note that rationalisability can easily be decided for a single square root of a one-variable polynomial (p = N sqrt = 1) based on the degree deg q 1 of the polynomial: the square root q 1 (x 1 ) can be rationalised (and the corresponding function ψ can be constructed explicitly) if and only if deg q 1 ≤ 2. The question of the rationalisability for p > 1 is much more involved, even in the presence of a single square root, see refs. [41][42][43][44][45].
In the following we discuss two special cases of eq. (3.4), in which we can express the Chen iterated integrals in terms of other classes of iterated integrals: • If N sqrt = 0, the Chen iterated integrals can be expressed in terms of MPLs evaluated at algebraic arguments.
• If N sqrt = 1 and deg q 1 = 3 or 4, the Chen iterated integrals can be expressed in terms of elliptic generalisations of MPLs evaluated at algebraic arguments.
In particular, we explain how we can algorithmically rewrite all Chen iterated integrals that meet these criteria in terms of MPLs and their elliptic analogues. This algorithm is well known, but we document it here because it will be an important tool to obtain compact analytic expressions for some master integrals that contribute to Bhabha scattering at two loops. Two comments are in order: First, in the presence of (multiple) square roots, it may be possible to change variables and rationalise (some of) the roots. In this way it may possible to reduce the problem to a situation covered by the criteria above, even though the original problem did not satisfy these conditions. Second, we stress that the aforementioned conditions are only sufficient, but by no means necessary, to rewrite Chen iterated integrals in terms of (elliptic) MPLs. In particular, there are several examples of Feynman integrals whose alphabets involve non-rationalisable square roots, but nevertheless, it was possible to express them in terms of ordinary MPLs, cf., e.g., refs. [24,46,47].

Rational alphabets without square roots
If the alphabet does not contain any square roots, or if all square roots can be rationalised, we can always express the Chen iterated integral in terms of MPLs evaluated at algebraic arguments. First, we can use the additivity of the logarithm to assume that all letters R i are irreducible polynomials. Next, we can use the homotopy-invariance of the integral to deform the contour γ into a new contour γ 0 with the same end-points x 0 and x 1 , without changing the value of the integral (except for picking up residues when we cross a pole).
We choose the new contour as follows: We order the integration variables in some way, which for simplicity we assume to be the natural order x 1 , . . . , x p . The contour γ 0 is then obtained as the concatenation of the straight-line segments γ r , 1 ≤ r ≤ p, defined by (3.5) Next we can iteratively apply the path-composition formula for iterated integrals, The previous equation allows us to reduce the integral to a linear combination of products of Chen iterated integrals over the straight-line segments γ r . On the line segment γ r the d log-forms take a particularly simple form. Indeed, since all the letters R i are polynomials in x, it is easy to see that R i (ϕ r (t)) is a polynomial in t. In the following we write where d denotes the degree of the polynomial R i (ϕ r (t)) and the quantities c j , 0 ≤ j ≤ d are algebraic functions of x 0 and x 1 . The pull-back of the one-form d log R i (x) to the straight-line segment γ r then reads The previous equation makes it manifest that on each straight-line segment the Chen iterated integral evaluates to ordinary MPLs with algebraic arguments c j , cf. eq. (3.2).

Alphabets with a single elliptic square root
We have already seen that a square root of a polynomial of degree three or four cannot be rationalised. More precisely, consider the set of points (x, y) constrained by where P n (x) is a polynomial of degree n in x. If n = 3 or 4, this equation defines an algebraic variety called an elliptic curve. It is thus not surprising that the case N sqrt = 1 and q 1 = P n , with n = 3 or 4, leads to generalisations of MPLs related to elliptic curves. We start by giving a lightning review of elliptic multiple polylogarithms (eMPLs), before we comment on the generalisation of the algorithm from section 3.1. More details can be found in appendix B. We focus here on the case n = 4, and we assume that P 4 has the form . Elliptic multiple polylogarithms (eMPLs) can then be defined as the iterated integrals [25,[48][49][50] with n i ∈ Z and c i ∈ C ∪ {∞}, and a = (a 1 , a 2 , a 3 , a 4 ) is the vector of the four branch points a i . Just like ordinary MPLs, eMPLs have at most logarithmic singularities, but no poles. The number n = k i=1 |n i | is called the weight of the eMPL, and the number of integrations k is its length. In the case where all the indices A i = ( n i c i ) are equal to ±1 0 , the integral is divergent and requires a special treatment similar to the case a n = 0 for ordinary MPLs, cf. eq. (3.3). We refer to appendix B for details. Also note that we need to be careful about how we choose the branches of the square root y = P 4 (x). We will come back to this point in section 5.
There are infinitely many integration kernels Ψ n (c, x, a) for given (c, x, a) in eq. (3.11). In concrete applications only a finite number of these kernels appear. Here we only list the kernels with |n| ≤ 1, which are of direct relevance to this paper. For n = 0, we have with c 4 = 1 2 √ a 13 a 24 , a ij = a i − a j . ω 1 is one of the two periods associated to the elliptic curve defined by the equation y 2 = P 4 (x), where K is the complete elliptic integral of the first kind (3.14) For |n| = 1, we have (with c = ∞) where we introduced the shorthand y c = P 4 (c). These functions have the property that the differential form dx Ψ ±1 (c, x, a) has a simple pole at x = c, and no other poles (except for dx Ψ 1 (c, x, a), which always has a pole at x = ∞). Note that the kernel Ψ 1 is identical to the kernel that defines ordinary MPLs (cf. eq. (3.2)), and so ordinary MPLs are a subset of eMPLs, The functions Z 4 (c, a) and G * ( a) in eq. (3.15) are in general transcendental functions. They can be expressed in terms of incomplete elliptic integrals of the first and second kind (see appendix B). Depending on the value of the argument c, in many applications it is possible to evaluate Z 4 (c, a) and G * ( a) in the form A + B iπ ω 1 , where A and B are algebraic functions of a and c [50]. In particular, in all cases relevant for this paper, Z 4 (c, a) and G * ( a) will always be algebraic functions (see section 5.1).
Let us now assume that the alphabet is given by d log-forms with arguments where r i,0 (x) are polynomials in x and q 1 (x) is a non-constant polynomial of degree at most four. We assume that the r i,j (x) do not have any common zero (otherwise we could factor out a term from this letter in the alphabet), and q 1 (x) does not have a double zero (otherwise we could factor it out of the square root). We allow r i,1 (x) to be zero (in which case this letter does not depend on the square root), but r i,0 (x) is assumed not to vanish. We can separate all the letters into even and odd parts: with Next we deform the integration path γ to the sequence of line segments γ r defined in eq. (3.5). Even letters give rise to integration kernels of the form dt t−c j and can be dealt with exactly as in section 3.1. We will not discuss them any further. For the odd letters, note that on the straight-line segment γ r , q 1 (ϕ r (t)) is a polynomial of degree n in t, with n ≤ 4. We write: a r = (a r,1 , . . . , a r,n ) . If n ≤ 2, then the square root P (r) (t) can be rationalised, and all letters are rational on the line segment γ r , and we do not need to discuss this case anymore. If the degree of P (r) (t) is three or four, then we cannot rationalise the square root on γ r . Instead, we obtain where the parity of R − i implies that F − r,i (t) is of the form where S r,i (t) is a rational function in t. At this point we make an important observation: By construction, the differential form d log R − i (ϕ r (t)) has only logarithmic singularities. This implies that S r,i (t) has poles of order at most one (possibly including a simple pole at infinity). As a consequence, S r,i (t) may be written as a linear combination of the functions Ψ −1 (c j , t, a r ) and Ψ 0 (c j , t, a r ). Hence, we can evaluate the integral on γ r in terms of eMPLs for the elliptic curve associated to P (r) (t). Note that, a priori, we may obtain eMPLs with different values of a r for each segment γ r .

Integration of the differential equations in terms of MPLs
After the general remarks of the previous section, we go back to the explicit solution for the system in eq. (2.5). The standard way to rationalize the first two square roots in eq. (2.6) is to turn to the dimensionless Landau variables x and y related to the Mandelstam invariants by −s In terms of these variables the Euclidean region s, t < 0 corresponds to 0 ≤ x, y ≤ 1 (using the symmetry of eq. (4.1) under (x, y) → (1/x, 1/y)). There exist different strategies to solve the differential equation in eq. (2.5). The first method consists in evaluating numerically the -expansion of the path-ordered exponential in eq. (2.9). This can be achieved by application of the Frobenius method to look for powerseries solutions of ordinary differential equations in the vicinity of regular singular points. The Frobenius method has been successfully used in various one-dimensional problems in the past [51][52][53][54], and more recently a generalisation of this method has been proposed in ref. [55] to deal with complicated multidimensional problems, see for example [56,57]. This strategy has also been implemented into the public code code DiffExp [58]. 3 The input data for this code are the matrices that define the differential equations, and the boundary conditions in some limit, e.g., at some point (x 0 , y 0 ). The code then uses this input to evolve the solution numerically from the point (x 0 , y 0 ) to some generic point (x, y). We fix the boundary conditions in the limit s, t → 0, which corresponds to (x 0 , y 0 ) = (1, 1) in terms of the Landau variables in eq. (4.1). Using the expansion by regions strategy implemented in the public code asy.m [59,60] (which is now included in the code FIESTA [61]) we obtain the following leading order asymptotic behaviour in this limit , for all the other elements. Let us emphasise that, in order to profit maximally from the automated code DiffExp, it is crucial that the input data are in an optimal form. This includes providing differential equations in canonical form and fixing the boundary conditions in a simple point (s = t = 0, i.e., x = y = 1). With our input, the code works very well and provides the possibility to obtain high-precision numerical results (100 digits accuracy and more), both in the Euclidean and the physical regions. In other words, DiffExp allows us to obtain high-precision numerical results for all master integrals and for all values of the input parameters. The code runs fast enough to allow its usage for practical applications. This is then of course sufficient for all phenomenological applications one has in mind. However, both from a formal and from a practical point of view, it may still be desirable in some situations to have full-fledged analytic representations of the master integrals in terms of functions that are well studied in the literature, e.g., MPLs. In the remainder of this section we explain how we can obtain analytic results for all master integrals but f 14 in terms of MPLs evaluated at algebraic arguments. The reason why f 14 is different will be discussed in section 5. Here it suffices to say that the square root r u only enters the differential equation for f 14 . Since our strategy of solving the differential equations in terms of MPLs will be closely based on the ideas from section 3.1, it is not suprising that an additional square root implies that the sufficient condition of section 3.1 is not satisfied. We will discuss the case of f 14 in detail in section 5, and we focus for the rest of this section only on the other master integrals.
In order to solve the master integrals in terms of MPLs, we start by observing that the square root r st does not appear when solving the differential equations up to weight 3 for all elements but f 37 , and at weight 4 for all elements but f i , i ∈ {35, 36, 37, 38, 39, 41, 43}.
The equations can be solved by integrating first in x, resulting in MPLs of the form G( c; x) with c i ∈ {0, −1, 1, −y, −1/y}. This allows us to fix the solution up to a function of y, which can be determined by substituting the results of the x-integration into the differential equations in y, checking that the variable x disappears in them, and finally solving them in terms of MPLs of the form G( c; y) with c i ∈ {0, −1, 1}, i.e., harmonic polylogarithms [37]. At this point the solution is fixed up to a set of undetermined integration constants, which we fix using the boundary conditions in eq. (4.2).
To evaluate f 37 at weights 3 and 4 and f i with i ∈ {35, 36, 38, 39, 41, 43} at weight 4, we have to deal with the square root r st . It can be rationalized by the following further change of variables 3) The resulting system of differential equations for 42 elements can be solved, first in w and then in y. Solving in w gives a linear combination of MPLs G(. . . ; w) with the letters the alphabet {l w 1 , . . . , l w 15 }, where l w i , i = 1, . . . , 11, are taken from the set 1, −1, − 2w 2 y 2 − y + 1 2 y 6 − 7y 5 + 12y 4 − 11y 3 + 12y 2 − 7y + 1 − 2w y 2 − y + 1 2 y 6 − 2y 5 + 4y 4 − 12y 3 + 4y 2 − 2y + 1 + 2y 10 − 11y 9 + 30y 8 − 62y 7 + 90y 6 − 89y 5 + 90y 4 − 62y 3 + 30y 2 − 11y + 2 .  Following this procedure, we obtain complete analytical results for f 37 at weights 3 and 4, and for the elements f i with i ∈ {35, 36, 38, 39, 41, 43} at weight 4, in terms of MPLs in w or y with the letters defined above. We have checked our analytic results for the master integrals against FIESTA [61] using also the forthcoming new release [62]. It turns out, however, (as it is easy to imagine from the dimension of the alphabet) that our results at weight 4 for these elements are rather complicated due to a very intricate branch-cut structure. In particular, evaluating them in this form at a phase-space point with GiNaC [63,64] meets problems connected both with timing and stability, so that our results for the complicated elements become impractical, and for phenomenological applications the numerical solution obtained using DiffExp might be preferable. For the same reason, we also recommend to turn to DiffExp also for the contribution of element 37 of weight 3. Of course, this does not imply that a better analytical representation in terms of MPLs does not exist, whose numerical evaluation could be much faster than automated numerical codes. Our analytic results obtained by direct integration of the differential equations could be used as a starting point to look for such a better analytic representation, if required.
Due to the symmetry f 14 (x, y) = f 14 (y, x), the differential equation with respect to y can easily be obtained by exchanging the roles of x and y in eq. (5.1). The square root in eq. (5.1) is identical to the square root that has appeared in ref. [16] in the computation of the first planar two-loop family for Bhabha scattering. The polynomial under the square root has degree four in both x and y, so it cannot be rationalised in one variable only.
However, this is not sufficient to exclude that one cannot rationalise it by performing a change of variable involving simultaneously x and y. In ref. [40] it was shown that the algebraic variety defined by ξ 2 = ∆(x, y) is a K3 surface. As a consequence, it cannot be rationalised by any rational change of variables, and the strategy of rationalising all square roots and evaluating all iterated integrals in terms of MPLs using ideas from section 3.1 cannot be applied. We can of course also obtain numerical results for f 14 by solving the canonical system with DiffExp (as described in the previous section), but it would of course be interesting to obtain analytic results also for f 14 .
We observe that the structure of the square root matches precisely the criteria of section 3.2. Indeed, the polynomial under the square root in eq. (5.1) has degree four in both x and y. So, if we keep one variable fixed, the square root defines an elliptic curve in the other. This reflects the fact that the K3 surface defined by the square root is elliptically fibered [40]. As a consequence, we can choose the integration path in eq. (3.5) and solve the differential equation in terms of eMPLs. Before we do this, we have to make an important comment. While the strategy of section 3.1 to solve eq. (5.1) does not apply, it does not exclude that a solution in terms of MPLs evaluated at algebraic arguments exists. On the contrary, we have found that it is possible to evaluate f 14 from its Feynman parameter representation using direct integration techniques (cf., e.g., refs. [65][66][67][68][69][70][71]). 4 The resulting expression, however, is extremely lengthy and involves MPLs evaluated at complicated algebraic arguments. We have not been able to confirm the final expression numerically, as its size and the complexity of the branch cuts renders the evaluation of the MPLs extremely challenging. Our representation obtained from direct integration is therefore not useful for practical purposes. As we will show now, the representation in terms of eMPLs is very compact.
Let us now explain how we can solve the system of differential equations forf in terms of MPLs and eMPLs. We first of all introduce new variables (x,ȳ) = (1 − x, 1 − y), and use PolyLogTools [71] to express all MPLs of the form G(. . . ; x) or G(. . . ; y) in eq. (5.1) in terms of G(. . . ;x) and G(. . . ;ȳ) respectively. For example, we find: We know from eq. (4. 2) that f 14 must vanish for s = t = 0, and sof (x = 0,ȳ = 0) = 0. Our strategy is then as follows. We first use the differential equation inȳ to evolvef from (x,ȳ) = (0, 0) to (x,ȳ) = (0,ȳ 0 ). We then use the differential equation inx to evolve from (x,ȳ) = (0,ȳ 0 ) to the generic point (x,ȳ) = (x 0 ,ȳ 0 ). In the following we assume 0 <x 0 ,ȳ 0 < 1 for concreteness (which corresponds to he Euclidean region s, t < 0). On the linex = 0 we find Hence, the square root disappears in the limitx → 0, and we can solve the differential equation on the linex in terms of MPLs. We find (relabellingȳ 0 asȳ): f (x = 0,ȳ) = 2π 2 log 2 − 3ζ 3 G 1 (ȳ) − π 2 G 1,1 (ȳ) + 2π 2 G 1,2 (ȳ) Next, let us discuss the solution on the lineȳ =ȳ 0 . Clearly, ∆(x, y) is not a perfect square for x = 0. Looking at the constraint ξ 2 = ∆(x, y 0 ) as a function of x with y 0 fixed, we see that it defines an elliptic curve. The branch points, i.e., the zeroes inx of We can use eq. (3.16) to interpret every MPL of the form G( b(ȳ);x) as an eMPL of the form E 4 n c(ȳ) ;x, a(ȳ) . Moreover, we can express all the algebraic coefficients multiplying the MPLs in the differential equation in terms of the functions Ψ ±1 (c,x, a(ȳ)). For example, we find with ξ = ∆(1 −x, 1 −ȳ). In the process, we discover the relations: At this point we need to make an important comment about how we choose the branches of the square root ξ = 1 −x, 1 −ȳ). From eq. (5.6) we can see that for 0 < y < 3 − 2 √ 2, the four branch points are real and ordered according to a 1 (y) < a 2 (y) < a 3 (y) < a 4 (y). The branches of the square root are chosen according to 1 , a 2 (y) ≤ x < a 3 (y) , i , a 3 (y) ≤ x < a 4 (y) .
The resulting differential equation can be solved by quadrature, and all integrals over x can be performed in terms of eMPLs using eq. (3.11). The final result reads (we relabel again (x 0 , y 0 ) as (x, y)): 1 ;x, a − 12Li 4 (−y) − 12Li 4 (y) − 2Li 2 (y) log 2 y − 2Li 2 (−y) log 2 y + π 2 + 8Li 3 (−y) log y + 8Li 3 (y) log y − 2ζ 3 log y − We have checked numerically that our final result forf (x, y) is correct by comparing against Fiesta. The eMPLs in eq. (5.11) were evaluated numerically using an in-house code (a public numerical implementation of eMPLs into GiNaC exists [72]). Note thatf (x, y) is a pure function of uniform weight four [50]. We find it interesting that such a compact analytic expression in terms of eMPLs exists, while our MPL expression obtained by naive direct integration is prohibitively large.

A compact analytic result for the first planar family
Let us conclude this paper by applying the techniques from section 5.1 to the planar family in fig. 1 (a). In ref. [16] it was shown that all integrals in this family can be expressed in terms of 23 master integrals, and all but one master integral were evaluated in terms of MPLs. The remaining master integral that could not be expressed in terms of MPLs in ref. [16] is (see fig. 3): This integral was evaluated in terms of a rather lengthy combination of MPLs with algebraic arguments in ref. [24]. Here we show that, by using the same strategy as in section 5.1, we can obtain a very compact representation for B(s, t, m 2 ) in terms of the same type of eMPLs as in eq. (5.11). The integral in eq. (5.13) is finite in four dimensions. Let us define B(x, y) by where the Landau variables (x, y) have been defined in eq. (4.1). Our starting point is the differential equation satisfied by B(x, y) [16]: The functions g i appearing in the right-hand side of eq. (5.15) are pure combinations of MPLs of weight three, Q(x, y) = (x + y)(1 + xy) x 2 y + xy 2 − 4xy + x + y = ∆(x, y) x 2 y + xy 2 − 4xy + x + y , (5.18) and the initial condition is B(x = 1, y = 1) = 0. Following exactly the same steps as in section 5.1, we find the following very compact expression: 1 ;x, a + 64ζ 4 − 32ζ 2 Li 2 (y) + 16Li 4 (y) + 8ζ 2 log 2 y + 1 3 log 4 y .
We have again checked our result numerically. Note that again we find a pure function of uniform weight four.

Conclusions
In this paper we considered the computation of the second family of planar master integrals relevant for Bhabha scattering at two loops in QED including the full dependence on the electron mass. Our primary tool for their analytic study was the method of differential equations, augmented by the choice of a canonical basis. We described how we obtained our canonical basis and showed that four different square roots appear in the calculation of the 43 master integrals that comprise it. We also provided boundary conditions at s = t = 0 for all master integrals. Together with the matrices defining the differential equations, this input is sufficient to produce high-precision numerical results for all master integrals and for all kinematic regions using the public code DiffExp. We provide a Mathematica notebook that allows one to evaluate all master integrals numerically with DiffExp as ancillary material with the arXiv submission. We also considered the analytic solution of the differential equations. Interestingly, the four square roots never appear all at the same time, and the differential equations for all master integrals but one can be solved in terms of MPLs by rationalising three of the four square roots by suitable changes of variables. For the contributions up to weight three, this procedure leads to analytic results which we present with the arXiv submission. While conceptually straightforward, we find that this procedure generates rather involved analytical expressions for the weight four part of the master integrals. The analytic results are available in electronic form from the authors upon request.
In the last part of the paper, we focused on the analytic computation of the remaining master integral, whose canonical differential equations contain three different square roots which cannot be rationalised at the same time. An analytic calculation in terms of MPLs cannot be easily obtained by direct integration of the differential equations due to the non-rationalisable square root. Instead, we show that compact analytic expressions can be obtained algorithmically in terms of elliptic multiple polylogarithms. We applied this idea in detail to our problem and obtained in this way a very compact analytic expression for this integral. We also showed that a similar compact expression can be obtained for another planar integral relevant for two-loop Bhabha scattering, whose calculation in terms of (a lengthy combination of) MPLs had been considered some years ago.

B Elliptic multiple polylogarithms
In this appendix we collect some technical material related to eMPLs not reviewed in section 3.2.
In the case where all the indices A i = ( n i c i ) are equal to ±1 0 , the integral in eq. (3.11) is divergent and requires a special treatment and we define instead, The rather ad-hoc looking form of eq. (B.1) is determined essentially uniquely by requiring the regularised eMPLs to share the same algebraic and differential properties as their convergent analogues. We refer to ref. [50] for a detailed discussion. The functions Z 4 (x, a) and G * ( a) that appear in eq. (3.15) can be expressed in terms of incomplete elliptic integrals of the first and second kind [25,50]: ,