Supersymmetric Higgs boson pair production at hadron colliders

We study the pair production of neutral Higgs bosons through gluon fusion at hadron colliders in the framework of the minimal supersymmetric standard model. We present analytical expressions for the relevant amplitudes, including both quark and squark loop contributions, and allowing for mixing between the superpartners of left- and right-handed quarks. Squark loop contributions can increase the cross section for the production of two $\mathrm{CP}$-even Higgs bosons by more than two orders of magnitude, if the relevant trilinear soft breaking parameter is large and the mass of the lighter squark eigenstate is not too far above its current lower bound. In the region of large $\mathrm{tan}\ensuremath{\beta},$ neutral Higgs boson pair production might even be observable in the $4b$ final state during the next run of the Fermilab Tevatron collider.