Double Higgs boson production at $e^+e^-$ colliders in the two-Higgs-doublet model

We study the double Higgs boson production processes $e^+e^- \to hh f\bar{f}$ ($f\neq t$) with $h$ being the 125 GeV Higgs boson in the two-Higgs-doublet model with a softly-broken $Z_2$ symmetry. The cross section can be significantly enhanced, typically a few hundreds percent, as compared to the standard model prediction due to resonant effects of heavy neutral Higgs bosons, which becomes important in the case without the alignment limit. We find a strong correlation between the enhancement factor of the cross section and the scaling factor of the $hf\bar{f}$ couplings under constraints from perturbative unitarity, vacuum stability and current experimental data at the LHC as well as the electroweak precision data.