  In this paper we study the renormalization of the product of two operators
$O_1=-\frac{1}{4} G^{\mu \nu}G_{\mu \nu}$ 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.
