  A computation of inclusive cross sections for neutral Higgs boson production
through gluon fusion and bottom-quark annihilation is presented in the MSSM
with complex parameters. The predictions for the gluon-fusion process are based
on an explicit calculation of the leading-order cross section for arbitrary
complex parameters which is supplemented by higher-order corrections: massive
top- and bottom-quark contributions at NLO QCD, in the heavy top-quark
effective theory the top-quark contribution up to N$^3$LO QCD including a soft
expansion for the $\mathcal{CP}$-even component of the light Higgs boson. For
its $\mathcal{CP}$-odd component and the heavy Higgs bosons the contributions
are incorporated up to NNLO QCD. Two-loop electroweak effects are also
incorporated, and SUSY QCD corrections at NLO are interpolated from the MSSM
with real parameters. Finite wave function normalisation factors ensuring
correct on-shell properties of the external Higgs bosons are incorporated from
the code $\texttt{FeynHiggs}$. For the typical case of a strong admixture of
the two heavy Higgs bosons it is demonstrated that squark effects are strongly
dependent on the phases of the complex parameters. The remaining theoretical
uncertainties for cross sections are discussed. The results have been
implemented into an extension of the numerical code $\texttt{SusHi}$ called
$\texttt{SusHiMi}$.
