  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.
