  We study the transverse momentum ($q_T$) spectrum of the Higgs boson produced
via the annihilation of heavy quarks ($s,c,b$) in proton-proton collisions.
Using soft-collinear effective theory (SCET) and working in the five-flavour
scheme, we provide predictions at three-loop order in resummed perturbation
theory (N$^3$LL$'$). We match the resummed calculation to full fixed-order
results at next-to-next-to-leading order (NNLO), and introduce a decorrelation
method to enable a consistent matching to an approximate N$^3$LO (aN$^3$LO)
result. Since the $b$-quark initiated process exhibits large nonsingular
corrections, it requires special care in the matching procedure and estimation
of associated theoretical uncertainties, which we discuss in detail. Our
results constitute the most accurate predictions to date for these processes in
the small $q_T$ region and could be used to improve the determination of Higgs
Yukawa couplings from the shape of the measured Higgs $q_T$ spectrum.
