We calculate the three-loop Wilson coefficients of all physically relevant dimension-four operators, i.e. \( {G}_{\mu \nu}^a{G}^{a,\mu \nu} \) , \( {m}_i{\overline{q}}_j{q}_j \) and \( {m}_i{m}_j{m}_k^2 \) , in the short-distance expansion of the time-ordered product of a pair of gauge-singlet vector, axial-vector, scalar and pseudo-scalar currents. The results are given for a general non-Abelian gauge theory with arbitrary (compact semi-simple) gauge group and nf light fermion flavors (quarks) in a common arbitrary representation of the gauge group, which includes QCD as a special case. In particular, we allow for arbitrary flavor contents of each of the currents. For the axial-vector current the included contributions from so-called singlet diagrams are consistent with the one-loop axial anomaly.
