  We study the modification to di-Higgs production via gluon fusion within the
context of the Minimal Supersymmetric Standard Model(MSSM) and the
Next-to-Minimal Supersymmetric Standard Model(NMSSM) in the parameter space
allowed by current experimental and theoretical constraints, and also relevant
to the Large Hadron Collider(LHC) experiments in the near future. The
calculation is based on the analytical expression of the leading order Feynman
amplitudes (which includes both quark and squark loops). We separate the
di-Higgs production cross section into resonant, non-resonant, and interference
parts, in order to better understand the mechanisms that are responsible for
the modification to di-Higgs production rate in different regions of the
allowed parameter space. We also investigate the sensitivity of High-Luminosity
LHC (HL-LHC) to the di-Higgs production in these low energy supersymmetry(SUSY)
models. Furthermore, we examine the complementarity between di-Higgs searches
and direct searches for BSM particles and precision Higgs couplings
measurements at the HL-LHC. We found that the di-Higgs production cross section
can be enhanced significantly through resonant production. In the region where
the resonant production cross section is small, di-Higgs receives a moderate
enhancement due to the modifications in the Higgs couplings. In addition, there
is a strong correlation between di-Higgs production and single Higgs
production, and di-Higgs is not as sensitive as single Higgs at the HL-LHC.
