  We present accurate predictions for the inclusive production of a Higgs boson
in proton-proton collisions, via gluon-gluon fusion. Our calculation includes
next-to-next-to-leading order (NNLO) corrections in perturbative QCD, as well
as the resummation of threshold-enhanced contributions to
next-to-next-to-next-to-leading logarithmic (N$^3$LL) accuracy, with the
inclusion of the recently-determined three-loop constant coefficient (sometimes
referred to as N$^3$LL' accuracy).
  Our result correctly accounts for finite top, bottom and charm masses at
leading order (LO) and next-to-leading order (NLO), and includes top mass
dependence at NNLO. At the resummed level the dependence on top, bottom and
charm mass is accounted for at NLL, while only the top mass at NNLL. The
all-order calculation is improved by a suitable choice of the soft terms,
dictated by analyticity conditions and by the inclusion of subleading
corrections of collinear origin, which improve the accuracy of the resummation
away from the threshold region.
  We present results for different collider energies and we study perturbative
uncertainties by varying renormalization and factorization scales. We find
that, at current LHC energies, the resummation corrects the NNLO result by as
much as 20 % at $\mu_R=\mu_F=m_H$, while the correction is much smaller, 5.5 %,
at $\mu_R=\mu_F=m_H /\,2$. While the central value of NNLO+N$^3$LL result
depends very mildly on the scale choice, we argue that a more reliable estimate
of the theoretical uncertainty is found if the perturbative scales are
canonically varied about $m_H$.
