  We present the calculation of the cross section and invariant mass
distribution for Higgs boson pair production in gluon fusion at next-to-leading
order (NLO) in QCD. Top-quark masses are fully taken into account throughout
the calculation. The virtual two-loop amplitude has been generated using an
extension of the program GoSam supplemented with an interface to Reduze for the
integral reduction. The occurring integrals have been calculated numerically
using the program SecDec. Our results, including the full top-quark mass
dependence for the first time, allow us to assess the validity of various
approximations proposed in the literature, which we also recalculate. We find
substantial deviations between the NLO result and the different approximations,
which emphasizes the importance of including the full top-quark mass dependence
at NLO.
