  Precise information on the Higgs boson self-couplings provides the foundation
for unveiling the electroweak symmetry breaking mechanism. Due to the scarcity
of Higgs boson pair events at the LHC, only loose limits have been obtained.
This is based on the assumption that the cross section is a quadratic function
of the trilinear Higgs self-coupling in the $\kappa$ framework. However, if
higher-order corrections of virtual Higgs bosons are included, the function
form would dramatically change. In particular, new quartic and cubic power
dependence on the trilinear Higgs self-coupling would appear. To get this new
function form, we have performed a specialized renormalization procedure
suitable for tracking all the Higgs self-couplings in each calculation step.
Moreover, we introduce renormalization of the scaling parameter in the $\kappa$
framework to ensure the cancellation of all ultraviolet divergences. With the
new function forms of the cross sections in both the gluon-gluon fusion and
vector boson fusion channels, the upper limit of
$\kappa_{\lambda_3}=\lambda_{\rm 3H}/\lambda_{\rm 3H}^{\rm SM}$ by the ATLAS
(CMS) collaboration is reduced from 6.6 (6.49) to 5.4 (5.37). However, it is
still hard to extract a meaningful constraint on the quartic Higgs
self-coupling $\lambda_{\rm 4H}$ from Higgs boson pair production data. We also
present the invariant mass distributions of the Higgs boson pair at different
values of $\kappa_{\lambda}$, which could help to set optimal cuts in the
experimental analysis.
