Zero-bin subtraction and the $q_{\mathrm{T}}$ spectrum beyond leading power

In this paper, we present an algorithm to construct the qT distribution at NLO accuracy to arbitrary power precision, including the assembly of suitable zero-bin subtrahends, in a mathematically well-defined way for a generic choice of rapidity-divergence regularisation prescription. In its derivation, we divide the phase space into two sectors, the interior of the integration domain as well as the integration boundary, which we include here for the first time. To demonstrate the applicability and usefulness of our algorithm, we calculate the NNLP corrections for Higgs hadroproduction for the first time. We observe that our approximate NNLP-accurate qT spectra replicate the asymptotic behaviour of the full QCD calculation to a much better degree than the previously available results, both within the $q_{\mathrm{T}}\to 0$ limit and in the large-qT domain for all the involved partonic processes. While playing a minor role at larger transverse momenta, we show that the newly incorporated boundary contribution plays a vital role in the $q_{\mathrm{T}}\to 0$ limit, where any subleading power accuracy would be lost without them. In particular, our NNLP-accurate qT expansion can approximate the exact qT distribution up to $q_{\mathrm{T}}\sim30\,\text{GeV}$ at the percent level for rapidities $|Y_H|\lesssim 3$.