  In this paper we consider two-loop two-, three- and four-point diagrams with
elliptic structure in the case of two different masses $m$ and $M$. The latter
diagrams generally arise within NRQCD matching procedures and are relevant for
parapositronium decay and top pair production at threshold. We present the
obtained results in several different representations: series solution with
binomial coefficients, integral representation and representation in terms of
generalized hypergeometric functions. The results are valid up to
$\mathcal{O}(\varepsilon)$ terms in $d=4-2\varepsilon$ space-time dimensions.
In the limit of equal masses $m=M$ the obtained results are written in terms of
elliptic constants with explicit series representation.
