Two-loop master integrals for pseudo-scalar quarkonium and leptonium production and decay

We compute the master integrals relevant for the two-loop corrections to pseudo-scalar quarkonium and leptonium production and decay. We present both analytic and high-precision numerical results. The analytic expressions are given in terms of multiple polylogarithms (MPLs), elliptic multiple polylogarithms (eMPLs) and iterated integrals of Eisenstein series. As an application of our results, we obtain for the first time an analytic expression for the two-loop amplitude for para-positronium decay to two photons at two loops.


Introduction
Current approaches to the calculation of higher-order corrections in perturbative Quantum Field Theory (QFT) are based on the decomposition of an observable into an independent set of Feynman integrals, usually called master integrals. While the details of the decomposition strongly depend on the specific QFT, the master integrals only depend on the kinematics of the process and the loop order. The computation of these Feynman integrals is thus an interesting problem in its own right, providing one of the crucial ingredients for the calculation of higher-order corrections to various physical observables.
Over the last decades, a lot of effort has been put into developing new techniques aimed at the efficient evaluation of Feynman integrals, both analytical and numerical. This effort has been focused on developing a better understanding of the classes of functions that these integrals evaluate to. It can be shown very generically that Feynman integrals evaluate to iterated integrals [1], with a non-trivial branch-cut structure that is constrained by physical considerations [2]. These very generic arguments still leave a large space of functions to explore, but it is by now well known that many Feynman integrals evaluate to iterated integrals with a very specific underlying geometry. The simplest kind of functions one finds are multiple polylogarithms (MPLs) [3], which, roughly speaking, correspond to iterated integrals over rational functions. The mathematical properties of MPLs are well understood. In particular, several computational tools have been developed to work with this class of functions, both for their analytic manipulation and their numerical evaluation. It is however known that starting at the two-loop order a new class of iterated integrals appears, where the integration kernels involve square roots of cubic or quartic polynomials which define an elliptic curve. If a single elliptic curve is present, the iterated integrals can be written in terms of elliptic multiple polylogarithms (eMPLs) [4][5][6]. Despite many recent developments, our understanding of the analytic structure of eMPLs is not nearly as mature as in the case of MPLs, and their numerical evaluation is still very challenging. Under certain conditions eMPLs can be expressed in terms of iterated integrals of Eisenstein series [7,8], for which more advanced numerical evaluation strategies are known [9]. Computing sets of master integrals involving elliptic curves is at the forefront of the problems that can currently be tackled.
In this paper we consider the complete set of master integrals appearing in the calculation of the two-loop amplitudes describing the production or decay of a pseudo-scalar bound state of a pair of massive fermions of the same flavour, be they quarkonium or leptonium bound states. For instance, the set of master integrals we consider is sufficient to compute the next-to-next-to-leading order (NNLO) QCD corrections to the hadro-production of a pseudo-scalar (cc) bound state (usually called η c ) at the LHC, or the NNLO QED corrections to the decay of an (e + e − ) bound state (usually called para-positronium) into two photons. These processes are of great phenomenological interest. For instance, quarkonium production offers interesting opportunities to study the interplay between the perturbative and non-perturbative regimes of QCD [10][11][12][13][14]. The master integrals contributing to the NNLO corrections to such processes involve both MPLs and eMPLs and, while some of them were already known in the literature [15][16][17][18][19][20][21][22][23][24][25][26][27][28], the complete set was not yet available in analytic form.
The purpose of this paper is two-fold: we present both analytic expressions for a complete set of two-loop master integrals contributing to the processes discussed above, and high-precision numerical evaluations of the integrals which can then be used for phenomenological studies. Regarding the analytic calculation, we obtain analytic expressions for all integrals by direct integration of their parametric representation. As already noted, we find that the integrals can be expressed in terms of MPLs, eMPLs and iterated integrals of Eisenstein series. The set of elliptic Feynman integrals involves two different elliptic curves: one elliptic curve belongs to the same family as the sunrise integral [29][30][31][32][33][34][35][36], while the second is an elliptic curve that appears in certain master integrals for tt production at hadron colliders [22,37]. Importantly, these two elliptic curves appear in independent sets of master integrals. Regarding the numerical evaluation of the master integrals, we present high-precision numerical results valid to 1000 digits. These numbers are obtained by numerically solving systems of differential equations within two slightly different approaches. More specifically, the results valid to 1000 digits are obtained with the auxiliary mass flow method [38][39][40] as implemented in AMFlow [40], and they are validated with results obtained with the generalised power series expansion [41,42] as implemented in diffexp [42]. These high-precision numerical results allow us to identify relations between the coefficients in the Laurent expansion of the master integrals in the dimensional regulator ǫ using the PSLQ algorithm [43]. These relations are important to obtain more compact analytic results, and it would certainly be interesting to understand how they can be generated more systematically. Our results can be found in a Mathematica-readable format at ref. [44].
The calculation of the various amplitudes which can be written in terms of the set of master integrals we consider in this paper will be discussed in detail in a companion paper [45]. We nevertheless present here analytic results for the two-loop QED corrections to the decay of true parapositronium, which were first obtained in numerical form more than 20 years ago in refs. [46,47], adding to the small but increasing number of physical quantities involving elliptic integrals which are known in analytic form [48][49][50][51][52].
The paper is structured as follows. In section 2 we discuss the set of master integrals we will consider in this paper. We present two types of relations beyond integration-by-parts identities which arise in degenerate kinematic configurations such as the ones corresponding to the amplitudes we consider. In section 3 we describe the analytic computation of the master integrals by direct integration, and we characterise the solutions in terms of MPLs, eMPLs and iterated integrals of Eisenstein series. In section 4 we summarise the main steps we followed for the numerical evaluation of the master integrals, along with the checks we did on our results. We also obtain relations between the elliptic master integrals using the high-precision numerical evaluations. Finally, in section 5 we discuss our results for the two-loop amplitude for the para-positronium decay to two photons, before we present our conclusion and outlook in section 6.

Kinematics and conventions
We consider the master integrals required to compute the two-loop perturbative corrections to the production or decay of a pseudo-scalar bound state of two massive fermions. This bound state can be a quarkonium bound state, in which case we consider higher-order QCD corrections, or a leptonium bound state, in which case we consider QED corrections. For concreteness, the discussion of this section focuses on the production of such a bound state, but it is clear that it also holds for its decay.
The perturbative corrections to the production of a bound state of massive fermions are systematically accounted for by considering the corrections to the short-distance process where Q and Q are fermions of mass m Q and the initial-state particles ab can be two gluons (gg), two photons (γγ) or a photon and a gluon (γg). We will consider the process in eq. (2.1) at leading order in an expansion in the relative velocity v of the QQ pair in the bound-state rest frame. This amounts to equating the heavy-fermion momenta p 1 and p 2 , and the kinematics effectively degenerate to those of a three-point process. More explicitly, assuming k 1 and k 2 incoming and p 1 and p 2 outgoing, we have The usual Mandelstam variables associated with four-point kinematics become Besides the perturbative corrections we consider here, there are also corrections related to higher orders in an expansion in v (see, e.g., ref. [53] and references therein). All these contributions must in general be taken into account for phenomenological predictions for the production or decay of quarkonium or leptonium bound states, but they fall outside of the scope of this paper. We will focus on the calculation of the master integrals that contribute to the two-loop amplitudes for the process in eq. (2.1). For quarkonium states, this will allow us to compute the two-loop amplitudes for the production and decay of the colour-singlet pseudo-scalar state 1 S [1] 0 in both the gg and γγ channels, or the production and decay of the pseudo-scalar colour-octet state 1 S [8] 0 in the gg and γg channels. For leptonium states, our set of integrals is sufficient to compute the two-loop QED corrections to para-positronium decaying into two photons, for which a numerical computation has already been performed in refs. [46,47], or the equivalent process for (true) muonium or tauonium. We will return to the decay of para-positronium in section 5. The calculation of the two-loop corrections to the production and decay of a pseudo-scalar bound state, both in colour-singlet and colour-octet, will be described in detail in a companion paper [45].
To determine the set of master integrals (MIs) that contribute to the process in eq. (2.1), we use standard techniques for the calculation of scattering amplitudes (we refer the reader to ref. [45] for more details). Let us nevertheless highlight here a consequence of the degenerate kinematics of eqs. (2.2) and (2.3). For generic four-point kinematics at two loops there are nine independent scalar products involving at least one loop momentum, but due to the degenerate kinematics of eqs. (2.2) and (2.3) only seven of them are independent. This fact must be taken into account when mapping all integrals into topologies, and we use the program Apart [54] to implement partialfraction relations systematically. Once all integrals have been sorted into topologies, we reduce them to MIs using integration-by-parts (IBP) relations [55,56]. We perform the IBP reductions with publicly available codes such as FIRE [57], LiteRed [58] and Kira [59].
Since there are 7 independent scalar products involving at least one loop momentum, all master integrals can be embedded in topologies involving at most 7 propagators, i.e., they can be written as where the D i denote inverse propagators and the a i take integer values. We refer the reader to appendix A for the explicit representation of each master integral in the form of eq. (2.4). We consider the integrals in dimensional regularisation in d = 4 − 2ǫ dimensions, and we normalise the integration measure as where γ E = −Γ ′ (1) is the Euler-Mascheroni constant. We note that, as a consequence of the degenerate kinematics in eq. (2.3), all master integrals are single-scale integrals whose explicit dependence on m 2 Q can be determined from dimensional analysis. Each integral can then be written as a Laurent series in ǫ, where dim(m I ) is half of the mass dimension of the integral m I , the F (k) I are constants, and we used the fact that two-loop master integrals have at most quadruple poles in ǫ. The goal of this paper is to determine these constants for each of the master integrals, up to the order in ǫ required to compute the two-loop amplitudes for the processes mentioned previously.

Partial-fraction and triangle relations
As we already noted, out of the 76 master integrals obtained after IBP reduction, several are available in the literature. To further reduce the number of integrals we need to compute, we discuss here a set of special identities beyond IBP or symmetry relations. The first kind of relations follows from partial-fraction relations due to the degenerate kinematics in eqs. (2.2) and (2.3). The second kind is obtained from a relation between three-point functions with a special mass configuration. We stress that we have not tried to find all identities that go beyond IBP relations, and it would undoubtedly be interesting to find and study such relations in a more systematic way.

Partial-fraction relations
The fact that partial-fraction relations are useful within the framework of quarkonium physics is well known (see, e.g., refs. [60,61]). For example, we have already mentioned that partial-fraction relations play an important role when sorting integrals into topologies. In addition, there are also partial-fraction relations that relate integrals from different topologies. We discuss some examples of such relations in this section.
In order to illustrate how partial-fraction relations arise, let us consider the following one-loop integral, with the propagators We then note that D 1 + D 3 = 2D 2 , from which it follows that   Diagrammatically, (2.10) In summary, this partial-fraction relation relates a four-point function to simpler three-point functions. A similar approach can be used to generate relations between two-loop integrals. The propagators in eq. (2.8) satisfy the relation which can diagrammatically be represented as where the red and blue dots connect to the rest of the diagram. It is possible to derive similar partial-fraction relations for other combinations of propagators. We can then use them to find relations among master integrals from different topologies. Using this type of identities, we found the following relations: (2.13) The last two relations have been derived by combining partial fraction relations with IBP relations and therefore involve the dimension d = 4 − 2ǫ and the mass scale m Q . The question naturally arises whether one could systematically incorporate these partial fraction relations at intermediate steps in the IBP reduction system to find all possible relations over the different topologies.

Triangle relations
Another special identity follows from a relation between certain triangle integrals. Consider a oneloop triangle integral with external legs p 1 , p 2 and k 1 = p 1 + p 2 , with k 2 1 = 0 and p 2 i = 0. We furthermore assume that two of the propagators have the same mass. Explicitly, we consider , (2.14) with Introducing Feynman parameters and after some manipulations, we find that (2.16) If we change variables according to x = 1/y, we get: Since this relation holds for arbitrary p 2 1 and p 2 2 , it can be used to relate multi-loop integrals having this triangle as a sub-diagram. In particular, it follows that = . (2.19) These integrals are related to m 10 and m 17 by IBP relations (see fig. 1), and so we can use this identity to relate m 10 to m 17 .

Analytic results for the master integrals
In the previous section we have defined our set of master integrals, which is composed of a total of 76 MIs. Taking into account results available in the literature [15][16][17][18][19][20][21][22][23][24][25][26][27][28] and the relations discussed in section 2.2, there are 38 master integrals that we have to consider. This can be for various reasons: some integrals are unknown, some are known in more general kinematic configurations and the limit cannot be smoothly taken, some have been computed but are expressed in different classes of functions, and some are known but not to the required order in the Laurent expansion around ǫ = 0. The integrals we must consider can be classified as follows: there are 16 four-point integrals, In this section, we review our approach to the analytic computation of the MIs listed above and discuss some of their analytic properties. Since we are only interested in their contributions to the NNLO corrections to the processes discussed in section 2 [45], we only focus on the terms in the Laurent expansion around ǫ = 0 that contribute at this perturbative order. In particular three MIs actually are absent from the NNLO corrections: m 5 , m 25 and m 38 contribute to the two-loop amplitudes only at O(ǫ). We nevertheless compute the leading order of m 5 and m 38 both analytically and numerically, while for m 25 we only present numerical results for its leading order. 1 We recall that for each master integral m I , our goal is to compute the numbers F (k) I in eq. (2.6). These numbers are special functions evaluated at specific numerical values. We will encounter three types of special functions, which are all particular instances of iterated integrals [1]: 2 1. Multiple polylogarithms (MPLs) [3], which are iteratively defined by G(a 1 , · · · , a n ; G(a 2 , · · · , a n ; t) , with G(; x) = 1.

Iterated integrals of Eisenstein series
with I(; τ ) = 1 and the f i are Eisenstein series of weight k i for some subgroup Γ ⊂ SL(2, Z).
They are special cases of modular forms of weight k i , i.e., holomorphic functions on the upper half-plane such that Presenting a detailed account of these functions and their properties would go beyond the scope of this paper, and we refer instead to the relevant literature (cf., e.g., refs. [6,9,[62][63][64][65] and references therein, for a discussion in a physics context).
Note that a given integral may be expressible in terms of more than one of these classes of iterated integrals. Depending on the representation chosen, the results can be more or less easy to manipulate and evaluate. Indeed, the understanding of these different types of iterated integrals is currently not on the same footing. For example, we know how to systematically simplify expressions involving MPLs and how to numerically evaluate them very efficiently (see, e.g., refs. [62,[66][67][68][69][70][71], and references therein). This is however not (yet) the case for eMPLs and iterated integrals of Eisenstein series. While first public codes for the numerical evaluation of eMPLs exist [72], these codes are not nearly as efficient as in the MPL case. In some cases, it is possible to write eMPLs as iterated integrals of Eisenstein series [63], and in this representation the numerical evaluation is much more efficient [9,31], allowing us to reach a precision comparable to what can be achieved for MPLs.

Direct integration
We have evaluated all the MIs for which no results were available in the literature via direct integration. Our starting point is the parametric representation of the scalar Feynman integrals in terms of Feynman parameters: where N is some normalisation factor, d = 4 − 2ǫ, Σ 0 = n i=1 x i , and U and F are the usual Symanzik polynomials. The representation in eq. (3.6) can be integrated order by order in ǫ using direct integration. 3 We distinguish the two cases: • MPL case: we are able to perform all the integrations in eq. (3.6) in terms of MPLs, and consequently the integral can be expressed in terms of MPLs evaluated at algebraic arguments; • eMPL case: we are not able to perform all the integrations in eq. (3.6) in terms of MPLs, because the integrand depends on the square root of a polynomial of the fourth order, which describes an elliptic curve. In this case we exploit the strategy outlined in refs. [37,75] to obtain an analytic expression in terms of eMPLs. (We refer to ref. [6] for how to relate iterated integrals involving square roots of quartic polynomials to the eMPLs in eq. (3.2)).
The MPL case is vastly discussed in the literature [67,[76][77][78][79][80]. In a nutshell, one starts from the Feynman-parameter representation in eq. (3.6) and notices that the Feynman-parameter integrals are projective integrals. This observation leads to what is usually known as the Cheng-Wu theorem in physics [81], which states that the integral is invariant under some transformations of the argument of the delta function, such as: Using the freedom afforded by this transformation, together with different integration orderings, one then looks for a way to integrate over all Feynman parameters so that each integration can be written in terms of MPLs. This can be done using public codes such as HyperInt [67], PolyLogTools [65] or MPL [80].
For the MIs which do not admit a representation in terms of MPLs, we use the strategy introduced in refs. [37,75]. More precisely, we proceed exactly as in the MPL case, but, due to the appearance of a square root defining an elliptic curve, we cannot perform all integrations in terms of MPLs. Instead, we find a transformation as in eq. (3.7) and an integration ordering such that we can perform all but the last integration in terms of MPLs. The last integration will depend on the final Feynman parameter x and the square root of a polynomial of degree 4 in x (denoted by y) describing the underlying elliptic curve. Schematically, the last integration will be of the form: where R(x, y) is a rational function in x and y, and J MPL (x, y) is a linear combination of MPLs with numerical coefficients, where we note that the square root y can also appear in the arguments of the MPLs. The next step is to rewrite the MPLs that appear in the last integration as eMPLs. This can be done systematically as described in ref. [37], by iteratively differentiating the MPLs and integrating them back in terms of eMPLs. For our MIs, these rather straightforward steps are complicated by the appearance of spurious square roots at either the second to last or the last integration. By spurious roots, we mean roots which are not related to the square root defining the elliptic curve. In the remaining part of this subsection we discuss how we overcame this issue (we note that spurious roots can also appear in the calculation of integrals involving only MPLs, and they can be dealt with in the same way; we focus here on the elliptic case since we found it to be by far the most challenging one). Let us first discuss the MIs for which spurious square roots appear in the second to last integration. This was the case for m 3 , m 5 , m 7 , m 13 , and m 14 , for which we could not find a transformation of the type of eq. (3.7) and an integration ordering that allowed us to rewrite all but the last Feynman-parameter integration in terms of MPLs. In this scenario, the next-to-last integration is of the form: where S is a rational function in x 1 , x 2 and the spurious square roots, and I MPL (x 1 , x 2 ) is a linear combination of MPLs whose arguments are rational functions in x 1 , x 2 and the spurious square roots. If there is a single spurious square root, there is a simple solution: we find a change of variables that rationalises the spurious root, for example with the approaches described in refs. [82][83][84], and then proceed with the integration over the new variables with the general strategy outlined above. 4 This was sufficient for MIs m 5 , m 7 , m 13 , and m 14 . Integral m 3 proved to be more challenging.
With the most convenient choice of transformation of the type of eq. (3.7) and most convenient integration order that we could find, we ended up with multiple spurious square roots appearing inside the same term. We found different terms involving different pairs of spurious roots, each still depending on two integration variables. In principle, one could attempt to rationalise all these square roots simultaneously, but this typically leads to very complicated changes of variables, which generate long expressions that are difficult to manipulate. Instead, we found that the expression could be simplified by inserting whereΣ is linear in the Feynman parameters and inx, into eq. (3.6): We then found that there was a choice of integration ordering that allows us to perform all but the last integral in terms of MPLs (more precisely, there were still spurious square roots, but all of them could be rationalised), bringing us into the general strategy discussed above.
To be more explicit, let us look more closely at the case of m 3 . Its Feynman parameter representation is where the Symanzik polynomials are: (3.13) We first note that by choosing Σ = x 1 +x 2 +x 4 we found the most compact intermediate expressions.
In the most convenient integration order that we could identify, we encounter the following spurious square roots where some of the spurious square roots appear simultaneously in the same terms. However, we found that by rewriting m 3 as with Σ as given above, and choosing the integration order we obtain an expression with two spurious square roots in the integration variables, which do not appear simultaneously within the same term. As such, we can proceed and rationalise the square roots in each term separately, and perform the next-to-last integration in terms of MPLs. 5 Let us now discuss how we deal with spurious square roots in the last integration step, which is a much more common occurrence. These spurious roots are generated when factorising the denominators in the second-to-last integration. In most cases, we find only spurious square roots involving a quadratic polynomial which can always be rationalised. However, for some integrals we encountered different spurious roots involving quartic or quintic polynomials, sometimes appearing in the arguments of the MPLs. This was for instance the case for m 3 , where we encountered more than 20 different spurious roots. In order to deal with these roots, we devised an algorithm that allows us to eliminate all spurious roots systematically. The key idea is that it is not necessary to bring individual terms into a form that can be integrated in terms of MPLs or eMPLs, but only the entire integrand. Our algorithm is iterative in the transcendental weight and proceeds as follows: 1. We first consider all terms of transcendental weight one. We then group them according to the x-dependence in their prefactor, where x denotes the last integration variable. Let us call g 1,i one of these combinations of MPLs of weight one with the x-dependent prefactor removed. We compute the derivative d dx g 1,i and construct the dlog-kernels. At this stage, all spurious roots in d dx g 1,i must disappear, otherwise they would not be spurious. We then integrate d dx g 1,i back in x and fix the boundary condition either analytically or with PSLQ. We repeat this procedure for all other weight one terms g 1,j .
2. We next consider all terms of transcendental weight two. The first step is again to group them according to the x-dependence of their prefactor. When computing the derivative of one such term, d dx g 2,i , we now have MPLs of weight one which might still involve spurious roots. These are dealt with as in the previous step. Next, we construct the dlog-kernels for the weight two contributions, and again all spurious roots disappear. We then integrate d dx g 2,i back in x and the boundary condition is fixed as at weight one. We repeat this procedure for all other weight two terms g 2,j .
3. We proceed with the same steps for weight three and then weight four, where each time the procedure requires dealing with all the lower weight terms that are generated when taking derivatives.
The crucial idea behind this algorithm is that, while individual terms may have spurious roots, the combination of terms that have the same weight and share the same integration kernel cannot, otherwise the square roots would not be spurious. Once all spurious roots have been removed, we can perform the last integral in terms of eMPLs for all integrals where spurious square roots appeared.
Following the steps described in this subsection we obtain fully analytic results for all the master integrals that we are interested in. We also obtained expressions for the integrals that had already been considered previously in the literature [15][16][17][18][19][20][21][22][23][24][25][26][27][28], and we provide results for the complete set of master integrals that contribute to the amplitudes mentioned in section 2. The analytic results for the MIs can be very lengthy, so we make them available in Mathematica-readable form at ref. [44]. The file analytics.m contains a replacement list of the form I , and the right-hand side of eq. (3.17) incorporates these relations. An independent set 6 of F  can be expressed entirely in terms of MPLs up to weight four evaluated at algebraic arguments. As already mentioned, the algebraic properties and the numerical evaluation of MPLs are well understood. In particular, using the public implementation in GiNaC [69] we obtained high-precision evaluations of these MIs to hundreds of digits. We then used the PSLQ algorithm [43] to fit our results to a basis of transcendental numbers. We observe that all of the MIs in eq. (3.18), except m 8 , m 38 , m 50 and m 51 , only involve MPLs evaluated at sixth roots of unity. This space of transcendental numbers has dimension 88, and a basis is known [85,86]. This observation allows us to obtain very compact expressions for these integrals. The four remaining master integrals, m 8 , m 38 , m 50 and m 51 , involve fourth roots of unity or algebraic numbers of the form a + b √ 2, where a and b are rational numbers. In this cases it is possible to construct a (conjectural) basis for these transcendental numbers using the results from section 2.5 of ref. [87].
The remaining 22 MIs not in eq. (3.18) cannot be expressed in terms of MPLs alone, but they involve eMPLs. 7 eMPLs depend on an elliptic curve, which is defined either by the square root of a quartic polynomial or, equivalently, by a value of τ as in eq. (3.2). We find two different elliptic curves in our computation, defined by: where K(λ) denotes the complete elliptic integral of the first kind: where m i and n i are integers and the ξ i are irrational numbers that are not rational multiples of τ (x) . Instead, they can be expressed in terms of the complete and incomplete elliptic integrals of the first kind evaluated at special arguments. The complete elliptic integral of the first kind was already defined in eq. (3.20), while its incomplete version is defined by All the ξ i can be written as rational linear combination of thez i listed in appendix B. If all arguments of an eMPL have the form aτ + b, with a, b ∈ Q (i.e., whenever ξ i = 0 in eq. (3.23)), the eMPL can be expressed in terms of iterated integrals of Eisenstein series, defined in eq. (3.4), in an algorithmic fashion [63]. Writing integrals in terms of iterated integrals of Eisenstein series presents the advantage that we know a basis for this type of iterated integrals and also efficient techniques for their numerical evaluation. We identified seven MIs for which their analytic expressions involve iterated Eisenstein series for two distinct subgroups Γ ⊆ SL • 20 MIs can be expressed in terms of eMPLs associated to the elliptic curve defined by τ (a) , and iterated integrals of Eisenstein series for Γ 1 (6).
• 2 MI can be expressed in terms of eMPLs associated to the elliptic curve defined by τ (b) , and iterated integrals of Eisenstein series for the congruence subgroup Γ 1 (4). We note once more that we only provide numerical results for m 25 since it is not strictly speaking needed for the processes we are concerned with. Nevertheless, given that it is coupled to m 24 , we know that it can be written in terms of the same set of functions.

Numerical results for the master integrals
In the previous section we discussed how we could obtain fully analytic results for the master integrals and we briefly reviewed the class of special functions that appear in the answer. We also noted that the numerical evaluation of MIs integrating to MPLs was a solved problem, and even used high-precision evaluations to write them in bases of transcendental numbers. For the numerical evaluation of eMPLs and iterated integrals of Eisenstein series, we follow the strategy presented in refs. [9,88,89]. In a nutshell, the integration kernels g (n) (z, τ ) and f j (τ ) in eqs. (3.2) and (3.4) are expanded in a so-called q-expansion (a Fourier expansion with q = e 2πiτ ). This series is truncated at a finite order, which, together with the convergence properties of the series representation, determines the precision of the numerical evaluation. The coefficients of this series expansion are MPLs, which can be numerically evaluated with standard tools. In the case of iterated integrals of Eisenstein series, it is known how to transform the integrals into a representation where the q-expansion converges very rapidly, and we can evaluate the integrals to hundreds of digits in an acceptable amount of time. For eMPLs, however, the convergence is in general very slow which in practice means we can only get a very small number of digits. For all integrals that involve eMPLs we obtained numerical results valid to O(10) digits. We have validated all our results by comparing them to completely independent numerical evaluations obtained with pySecDec [90], with which one can achieve a similar level of precision.
While the efficiency of the numerical evaluation of eMPLs will surely improve in the future, alternative approaches can already be used to obtain high-precision numerical evaluations for MIs involving eMPLs. These evaluations can then be used in phenomenological applications. In this section we explain how we bypassed the evaluation of eMPLs to obtain numerical results for all master integrals with a precision of 1000 digits. This precision is more than sufficient for any phenomenological application, and also allows us to use the PSLQ algorithm to simplify and find relations between the various F (k) I , as will be discussed in section 4.2. We discuss two independent calculations, which have in common the fact that they are based on numerically solving systems of differential equations.

High-precision numerical results from differential equations
To obtain high-precision numerical evaluations for the elliptic MIs that involve eMPLs, we construct systems of differential equations for them [15,[91][92][93], which are then solved numerically. These calculations are performed with two alternative approaches. In the first one, we build the differential equations and provide the required boundary conditions ourselves. The system is then solved in terms of generalised power series [41,42] using diffexp [42]. In the second approach, we use AMFlow [40], which automatises the construction and solution of the system of differential equations [38][39][40]. The integrals involving eMPLs can be collected in four topologies (topologies 3, 4, 5 and 6 in the notation of appendix A). In topology 6, however, the only elliptic integrals are m 5 , which only contributes at order ǫ to the two-loop amplitudes of ref. [45], and m 15 , which can be written in terms of iterated integrals of Eisenstein series, and can thus be evaluated efficiently to hundreds of digits from its analytic representation. We have thus only evaluated the elliptic integrals of topology 6 with AMFlow, since the evaluation with diffexp is only used as a check in the other topologies.
Let us briefly discuss the calculation with diffexp, since AMFlow implements similar steps in an automated way. In a nutshell, if I is a basis of master integrals for a given topology, we compute the system of differential equations where x denotes all the kinematic variables on which the integrals I depend. We then numerically solve this differential equation order by order in ǫ. We note, however, that we cannot straightfor- wardly apply the differential-equation approach to the master integrals in figs. 1 to 4. Indeed, as noted in eq. (2.6), they depend on a single variable and as such they satisfy a trivial differential equation. Instead, we will consider two-scale versions of the master integrals and evaluate them at the point corresponding to the kinematics of eq. (2.3). We are assisted in this task by the fact that a generalisation of topologies 3 and 4 was considered in refs. [21,23] and a generalisation of topology 5 in ref. [22], precisely from the perspective of their differential equations. For each topology we setup the calculation in the following way: • For topology 3: we consider the topology defined by the diagram in fig. 5a, evaluated at λ 2 = 0. Our basis is similar to the one chosen in refs. [21,23]. 10 For completeness, all the information required to completely define the basis we chose can be found in the folder /diffexp/ of ref. [44].
• For topology 4: we consider the topology defined by the diagram in fig. 5b, evaluated at λ 2 = 0. The same comments as for topology 3 apply. The complete information to define our basis can be found in the folder /diffexp/ of ref. [44].
• For topology 5: we consider the topology defined by the diagram in fig. 5c, evaluated at λ 2 = 4m 2 Q . The basis we choose is that of ref. [22], up to some trivial normalisation factors. The complete information to define our basis can be found in the folder /diffexp/ of ref. [44].
Constructing the differential equations for these three sets of master integrals is straightforward with standard approaches [15]. In all cases, our bases are such that the differential equation matrices M ( x; ǫ) in eq. (4.1) are polynomials in ǫ (with our choice of bases, they are polynomials of degree 2). In the sub-sectors that only couple integrals which evaluate to MPLs, the differential equation is in canonical (dlog) form, that is, the corresponding entries of M ( x; ǫ) consist of dlog-forms and are proportional to ǫ. To solve the differential equations we use the initial conditions given in refs. [21][22][23] within diffexp. This allows us to numerically solve the three systems of differential equations to hundreds of digits at the relevant kinematic points (if the integrals are logarithmically divergent at this point, we instead compute a generalised series expansion). We then use IBP relations to relate the bases chosen to solve the differential equations to the master integrals in figs. 1 to 4.
Up to some differences in the way the boundary conditions are determined, the evaluation within AMFlow is based on the same procedure. The different steps are automated so that one can simply require the evaluation of the integral at a phase-space point. Using AMFlow, we were able to evaluate the integrals to very high precision (O(1500) digits), but to be conservative we keep only the first 1000 digits. We found that we could achieve a higher precision with AMFlow than with diffexp, so we quote the former as our final results. While the difference is irrelevant for phenomenological studies, the highest precision of the AMFlow evaluations were useful in identifying extra relations between different integrals using PSLQ.
We find complete agreement between the numerical evaluations within the two approaches, and also with the pySecDec evaluations. Furthermore, in this process we also reevaluate master integrals that do not involve eMPLs, which we find to completely agree with the high-precision numerical evaluation of their analytic expressions. Our results, correct to 1000 digits, can be found in the file numerics.m at ref. [44].

Relations among elliptic master integrals
The fact that the algebraic properties and the numerical evaluation of MPLs is well understood allowed us to use the PSLQ algorithm to express the MIs that evaluate to MPLs in a basis of MPLs with specific arguments. This results in rather compact analytic representations for those MIs. This is in contrast to the case of MIs that involve eMPLs and/or iterated integrals of Eisenstein series where, as already noted in section 3.1, analytic expressions can be extremely long. Indeed, it is currently very complicated to simplify expressions involving eMPLs, because not much is known about how to manipulate such functions.
In order to simplify our analytic expressions, we were nevertheless able to find relations between the coefficients in the Laurent expansion of certain elliptic master integrals. We were motivated by the fact that the poles of the two-loop amplitudes for quarkonium or leptonium production or decay considered in ref. [45] are free of elliptic contributions. 11 Using the PSLQ algorithm we found the following relations: Re G 0, 0, e −2iπ/3 , 1; 1 , Re G 0, 0, e −iπ/3 , −1; 1 + 99 320 Re G 0, 0, e −2iπ/3 , 1; 1 + 3Re G e −iπ/3 , 1, 1, −1; 1 + 3 2 Re G e −iπ/3 , 1, −1; 1 log 2 , Re G e −iπ/3 , 1, −1; 1 log 2 ,  where ζ n denotes the Riemann zeta function evaluated at n (in particular, ζ 2 = π 2 /6 and ζ 4 = π 4 /90) and Cl 2 (x) denotes the Clausen function, Cl 2 (x) = Im Li 2 (e ix ) . It would be very interesting to find a systematic way to derive such identities, and some might follow from the type of relations discussed in section 2.2. In practice, using these relations we are able to write some of the lengthiest elliptic integrals in terms of much more compact expressions.
5 Analytic results for the para-positronium decay up to NNLO As mentioned in section 2, the master integrals we have discussed in this paper allow us to compute the NNLO contributions to the production and decay of several quarkonium and leptonium pseudo-scalar bound states. While we leave a more extensive discussion of such calculations to a separate publication [45], we finish this paper by discussing the NNLO corrections to the parapositronium (e + e − ) decay to two photons, presenting for the first time complete analytic results for this contribution. The para-positronium state is a pseudo-scalar particle with spectroscopic notation 0 +− (J P C ). 12 The most precise experimental measurement for its decay width includes the decay to any even number of photons, giving [95] Γ exp.
p-Ps decay = (7990.9 ± 1.7) (µs) −1 . (5.1) The leading-order of the decay to four photons is of the same order in perturbation theory as the NNLO corrections to the decay to two photons. We write where the leading-order contribution to Γ p-Ps→4γ has been computed analytically in ref. [96], The NNLO corrections to Γ p-Ps→γγ were first computed in purely numerical form more than 20 years ago [46,47] (see also ref. [97]), and in this section we present them for the first time in analytic form.
We can express the decay of the para-positronium to two photons up to NNLO accuracy in QED as [46,97,98], where α em is the electromagnetic fine structure constant and m e is the electron mass. The terms of the form α n em log k α em and K 2,soft are corrections related to the leptonium bound state [98][99][100][101]. These and the LO and NLO corrections were already known in analytic form in the literature. We note that K 2 and K 2,soft are independently divergent (each has a so-called Coulomb singularity) but their sum is finite. The master integrals that we have computed in this paper allow us to obtain for the first time the analytic results for the two-loop contribution K 2 .
We follow the presentation of refs. [46,47] and express the two-loop virtual amplitudes contributing to K 2 in terms of three different types of terms: (a) regular diagrams without fermion loops as shown in fig. 6a, (b) contributions with fermion loops in diagrams of type light-by-light scattering as in fig. 6b and (c) contributions with vacuum polarisations as in fig. 6c. These different contributions have ultraviolet poles which can be removed by renormalisation. After this process, a single pole remains which is related to a Coulomb singularity (see, e.g., ref. [46]). We then write: , denotes the two-loop amplitude obtained after subtracting the ultraviolet poles, which can be done contribution by contribution in the decomposition of fig. 6, and  where we have kept the elliptic contributions in the symbolic form F (k) I , whose analytic representation can be found at ref. [44]. We note that the two-loop amplitude exhibits functions of maximal weight four (for MPLs) and maximal length four (for elliptic functions), which is in full agreement with the conjectural property that scattering amplitudes at n-loops should exhibit functions of maximal weight/length 2n. Using our high-precision numerical evaluations, these different contributions evaluate to where the numbers are truncated to 20 digits after the decimal. The NNLO contribution K 2 in eq. (5.4) can then be decomposed as where K 1 = π 2 4 − 5 is the NLO contribution, and the remaining terms are twice the real part of eqs. (5.10), (5.11), (5.12) and (5.6) respectively. The soft contribution takes the form [46], Therefore, combining all these expressions, we obtain the high-precision numerical result where, as before, all number are truncated to 20 digits after the decimal. We note that the finite part of K 2 is dominated by the (negative) contribution of K 2,reg . The finite piece of K 2,soft has a similarly sized positive contribution, leading to a NNLO perturbative correction of O(1).
Since the precision of our numerical evaluation of K 2 depends on the numerical evaluation of our master integrals, which are correct up to over 1000-digits accuracy, they do not contribute to the theory error associated with the NNLO decay width. The theory error is determined by the precision to which the QED coupling α em and the electron mass m e are known. We use for the coupling α em = (7.2973525693 ± 0.0000000011) × 10 −3 , (5. 16) and for the electron mass m e = (0.5109989500 ± 0.00000000015) MeV, (5.17) where the errors are correlated with r = −0.99998 [102,103]. The propagation of these correlated errors to the decay width can be determined with: which is in full agreement with the experimental measurement quoted in eq. (5.1).

Conclusions
In this paper we have computed the complete set of two-loop master integrals contributing to the NNLO corrections to the production or decay of a pseudo-scalar bound state of two massive fermions of the same flavour. We have presented both analytic results as well as high precision numerical evaluations. The analytic results involve both MPLs and eMPLs (and iterated integrals of Eisenstein series). The presence of elliptic integrals presented the main challenge to the evaluation of this set of integrals. The high-precision numerical evaluations are comparatively simple to obtain with modern tools such as diffexp [42] and AMFlow [40]. Besides being important for phenomenological studies, these high-precision numerical evaluations also allowed us to considerably simplify the analytic expressions using the PSLQ algorithm. All our results can be obtained from the repository at ref. [44]. As an example of the type of processes these integrals can be used for, we recomputed the two-loop corrections to the decay of para-positronium into two photons [46,47,97], presenting both analytic results and high-precision numerical evaluations. Computing Feynman integrals involving elliptic functions is still a very challenging task. Moreover, once analytic expressions are obtained, it is not yet known how to simplify them or even how to efficiently evaluate them. We expect this situation to improve substantially in the coming years as more and more processes involving elliptic functions are computed, and the analytic results we obtained in this paper will be very useful in this process. For instance, understanding how to systematically find the relations in section 4.2 will be of great importance: in our calculation, it allowed us to rewrite the lengthiest elliptic integrals in terms of simpler ones. We also expect that there will be important developments in the numerical evaluation of eMPLs, and reproducing the high-precision numerical evaluations of the master integrals from their analytic representation will provide an important test.
The set of master integrals we have computed will open the door to new NNLO predictions for processes involving bound states of heavy fermions. Besides the para-positronium decay to two photons discussed in this paper, for which experimental precision is expected to increase at future experiments, see e.g. [104,105], and the equivalent process for (true) para-muonium and (true) paratauonium, our results also allow us to study the production and decay of quarkonium states at an unprecedented level of precision. An interesting prospect is to use charmonium production to study the gluon parton distribution function of the proton, as these are rather unconstrained at scales close to the mass of the charm quark. Furthermore, it is interesting to study the convergence of the perturbative expansion as the strong coupling α s is not so small at these scales, see refs. [106,107] and references therein. We will discuss several of these processes in a companion paper [45].
Regarding the evaluation of two-loop master integrals for processes with quarkonium or leptonium states, the next steps are clear. The set of integrals we considered here correspond to very simple kinematics, only depending on the mass of the heavy fermion. Adding an extra particle in the final state would lead to richer kinematic configurations and allow us to study other processes, for instance, hadro-production or photo-production of vector bound states with spectroscopic notation 3 S 1 , commonly called J/ψ (cc) and Υ (bb). Due to the Landau-Yang theorem, the LO of these production processes involves an additional gluon in the final state making it effectively a 2-to-2 process (gg → J/ψ + g). Similarly, we could study the p T -distribution of η Q at NNLO accuracy by computing the two-loop corrections to the process gg → η Q + g. Tackling such calculations will pose many challenges. In particular, the numerical evaluation will be much more involved as the integrals will depend on several variables. As such, a single high-precision evaluation of master integrals will no longer be sufficient.

B Arguments of the eMPLs
In this appendix we present the analytic expressions for the arguments of the eMPLs as defined in eq. (3.23). We find that we need six different values ofz i for the eMPLs associated to the elliptic curve τ (a) , and just one, which we denotez b , for the eMPLs described by τ (b) . Furthermore, we checked using PSLQ that there is no rational linear combination of the form where F(φ|m) was defined in eq. (3.24) and we set φ i = sin −1 α i , with (−1 + 2i) + √ 1 − 2i , Similarly, the numberz b is defined by the elliptic integral: We therefore havẽ where in this case φ b is defined as φ b = sin −1 α b , with