  The virtual two-loop corrections for Higgs production in gluon fusion are
calculated analytically in QCD for arbitrary Higgs and quark masses. Both
scalar and pseudo-scalar Higgs bosons are considered. The results are obtained
by expanding the known one-dimensional integral representation in terms of
m_H/m_q, and matching it with a suitably chosen ansatz of Harmonic
Polylogarithms. This ansatz is motivated by the known analytic result for the
Higgs decay rate into two photons. The method also allows us to check this
result and to extend it to the pseudo-scalar decay rate.
