  We consider the current experimental constraints on the parameter space of
the MSSM and NMSSM. Then in the allowed parameter space we examine the Higgs
pair production at the 14 TeV LHC via $b\bar{b}\to hh$ ($h$ is the 125 GeV
SM-like Higg boson) with one-loop SUSY QCD correction and compare it with the
production via $gg\to hh$. We obtain the following observations: (i) For the
MSSM the production rate of $b\bar{b} \to hh$ can reach 50 fb and thus can be
competitive with $gg \to hh$, while for the NMSSM $b\bar{b} \to hh$ has a much
smaller rate than $gg \to hh$ due to the suppression of the $hb\bar{b}$
coupling; (ii) The SUSY-QCD correction to $b\bar{b} \to hh$ is sizable, which
can reach $45\%$ for the MSSM and $15\%$ for the NMSSM within the $1\sigma$
region of the Higgs data; (iii) In the heavy SUSY limit (all soft mass
parameters become heavy), the SUSY effects decouple rather slowly from the
Higgs pair production (especially the $gg\to hh$ process), which, for $M_{\rm
SUSY}=5$ TeV and $m_A<1$ TeV, can enhance the production rate by a factor of
1.5 and 1.3 for the MSSM and NMSSM, respectively. So, the Higgs pair production
may be helpful for unraveling the effects of heavy SUSY.
