We calculate the unpolarized and polarized two–loop massless off–shell operator matrix elements in QCD to O(ε) in the dimensional parameter in an automated way. Here we use the method of arbitrary high Mellin moments and difference ring theory, based on integration-by-parts relations. This method also constitutes one way to compute the QCD anomalous dimensions. The presented higher order contributions to these operator matrix elements occur as building blocks in the corresponding higher order calculations up to four–loop order. All contributing quantities can be expressed in terms of harmonic sums in Mellin–N space or by harmonic polylogarithms in z–space. We also perform comparisons to the literature
