  We compute the next-to-leading virtual QCD corrections to the partonic cross
section of the production of a Higgs boson in association with a $Z$ boson in
gluon fusion. The calculation is based on the recently introduced method of
evaluating the amplitude via an expansion in terms of a small transverse
momentum. We generalize the method to the case of different masses in the final
state and of a process not symmetric in the forward-backward direction
exchange. Our analytic approach gives a very good approximation (better than
percent) of the partonic cross section in the center of mass energy region up
to $\sim 750 \,\textrm{GeV}$, where at the LHC $\sim 98\%$ of the total
hadronic cross section is concentrated.
