  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 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 $\tan \beta$, neutral
Higgs boson pair production might even be observable in the $4 b$ final state
during the next run of the Tevatron collider.
