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.