  We report a new calculation of the soft-gluon threshold resummation for the
Higgs boson pair production in the dominant production mode -- gluon-gluon
fusion -- up to the next-to-next-to-next-to-leading logarithmic (N$^3$LL)
accuracy. After matching N$^3$LL to the next-to-next-to-next-to-leading order
(N$^3$LO) QCD calculation in the infinite top quark mass approximation, we show
that the central values of the inclusive cross sections are quite stable with
respect to N$^3$LO, while the conventional renormalisation and factorisation
scale uncertainties are reduced by a factor of two, reaching to the subpercent
level. Our study further consolidates the good asymptotic perturbative
convergence. After combining with the full top-quark mass dependent
next-to-leading order QCD results, our most advanced predictions are presented
for both the inclusive total cross sections and the differential invariant mass
distributions of the Higgs pair.
