In this paper we study the renormalization of the product of two operators $O_1=-\frac{1}{4} G^{����}G_{����}$ in QCD. An insertion of two such operators $O_1(x)O_1(0)$ into a Greens function produces divergent contact terms for $x\rightarrow 0$. In the course of the computation of the operator product expansion (OPE) of the correlator of two such operators $i\int\!\mathrm{d}^4x\,e^{iqx} T\{\,O_1(x)O_1(0)\}$ to three-loop order we discovered that divergent contact terms remain not only in the leading Wilson coefficient $C_0$, which is just the VEV of the correlator, but also in the Wilson coefficient $C_1$ in front of $O_1$. As this correlator plays an important role for example in QCD sum rules a full understanding of its renormalization is desireable. This work explains how the divergences encountered in higher orders of an OPE of this correlator should be absorbed in counterterms and derives an additive renormalization constant for $C_1$ from first principles and to all orders in perturnbation theory. The method to derive the renormalization of this operator product is an extension of the ideas of a paper by Spiridonov and can be generalized to other cases.

v2: this is the version accepted by JHEP; more detailed discussion of phenomenological applications