Two-loop QED Operator Matrix Elements with Massive External Fermion Lines

The two-loop massive operator matrix elements for the fermionic local twist--2 operators with external massive fermion lines in Quantum Electrodynamics (QED) are calculated up to the constant terms in the dimensional parameter $\epsilon = D - 4$. We investigate the hypothesis of Ref. \cite{BBN} that the 2--loop QED initial state corrections to $e^+e^-$ annihilation into a virtual neutral gauge boson, except power corrections of $O((m_f^2/s)^k), k \geq 1$, can be represented in terms of these matrix elements and the massless 2-loop Wilson coefficients of the Drell-Yan process.