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.