  We compute the two-loop massless QCD corrections to the four-point amplitude
$g+g \rightarrow H+H$ resulting from effective operator insertions that
describe the interaction of a Higgs boson with gluons in the infinite top quark
mass limit. This amplitude is an essential ingredient to the third-order QCD
corrections to Higgs boson pair production. We have implemented our results in
a numerical code that can be used for further phenomenological studies.
