  Understanding the Higgs potential by measuring its self-interactions is
fundamental in answering several big questions, such as electroweak symmetry
breaking, electroweak baryogenesis, electroweak phase transition, and
electroweak vacuum stability. The most promising way to probe the Higgs
potential is to detect Higgs boson pair final state at high-energy colliders.
In this talk, we report a recent perturbative calculation for the di-Higgs
gluon-fusion process by taking into account N$^3$LO QCD radiative corrections
in the approximation of infinite top quark mass limit. Finite top quark mass
effects are also incorporated with several approximate schemes, which are known
to be crucial in phenomenological applications. We show a very good asymptotic
perturbative convergence at $\mathcal{O}(\alpha_s^5)$, and demonstrate that the
remaining scale uncertainty is only at percent level.
