  The Higgs boson pair production via gluon fusion at high-energy hadron
colliders, such as the LHC, is vital in deciphering the Higgs potential and in
pinning down the electroweak symmetry breaking mechanism. We carry out the
next-to-next-to-next-to-leading order (N$^3$LO) QCD calculations in the
infinite top-quark mass limit and present predictions for both the inclusive
and differential cross sections, albeit the differential distributions other
than the invariant mass distribution of the Higgs boson pair are approximated
at N$^3$LO. Such corrections are indispensable in stabilising the perturbative
expansion of the cross section in the strong coupling $\alpha_s$. At the
inclusive level, the scale uncertainties are reduced by a factor of four
compared with the next-to-next-to-leading order (NNLO) results. Given that the
inclusion of the top-quark mass effects is essential for the phenomenological
applications, we use several schemes to incorporate the N$^3$LO results in the
infinite top-quark mass limit and the next-to-leading order (NLO) results with
full top-quark mass dependence, and present theoretical predictions for the
(differential) cross sections in the proton-proton collisions at the
centre-of-mass energies $\sqrt{s}=13,14,27$ and $100$ TeV. Our results provide
one of the most precise theoretical inputs for the analyses of the Higgs boson
pair events.
