Two-loop massless QCD corrections to the $g+g\to H+H$ four-point amplitude

We compute the two-loop massless QCD corrections to the four-point amplitude $g+g\to 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.