  An important ingredient for the calculation of Higgs boson properties in the
infinite top quark mass limit is the knowledge of the effective coupling
between the Higgs bosons and gluons, i.e. the Wilson coefficients $C_H$ and
$C_{HH}$ for one and two Higgs bosons, respectively. In this work we calculate
for the first time $C_{HH}$ to four loops in a direct, diagrammatic way,
discussing in detail all issues arising due to the renormalization of operator
products. Furthermore, we also calculate the Wilson coefficient $C_H$ for the
coupling of a single Higgs boson to gluons as well as all four loop decoupling
relations in QCD with general SU$(N_c)$ colour factors. The latter are related
to $C_H$ and $C_{HH}$ via low-energy theorems.
