  In this Letter, we present for the first time a calculation of the complete
next-to-leading order corrections to the $gg \to ZH$ process. We use the method
of small mass expansion to tackle the most challenging two-loop virtual
amplitude, in which the top quark mass dependence is retained throughout the
calculations. We show that our method provides reliable numeric results in all
kinematic regions, and present phenomenological predictions for the total and
differential cross sections at the Large Hadron Collider and its future
upgrades. Our results are necessary ingredients towards reducing the
theoretical uncertainties of the $pp \to ZH$ cross sections down to the
percent-level, and provide important theoretical inputs for future precision
experimental collider programs.
