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 ε=D−4. We investigate the hypothesis of Berends et al. (1988) [1] that the 2-loop QED initial state corrections to e+e− annihilation into a virtual neutral gauge boson, except power corrections of O((mf2/s)k), k⩾1, can be represented in terms of these matrix elements and the massless 2-loop Wilson coefficients of the Drell-Yan process.