  Higgs-pair production via gluon fusion is the dominant production mechanism
of Higgs-boson pairs at hadron colliders. In this work, we present details of
our numerical determination of the full next-to-leading-order (NLO) QCD
corrections to the leading top-quark loops. Since gluon fusion is a
loop-induced process at leading order, the NLO calculation requires the
calculation of massive two-loop diagrams with up to four different mass/energy
scales involved. With the current methods, this can only be done numerically,
if no approximations are used. We discuss the setup and details of our
numerical integration. This will be followed by a phenomenological analysis of
the NLO corrections and their impact on the total cross section and the
invariant Higgs-pair mass distribution. The last part of our work will be
devoted to the determination of the residual theoretical uncertainties with
special emphasis on the uncertainties originating from the scheme and scale
dependence of the (virtual) top mass. The impact of the trilinear
Higgs-coupling variation on the total cross section will be discussed.
