# Fiducial distributions in Higgs and Drell-Yan production at N³LL+NNLO

Wojciech Bizoń, Xuan Chen, Aude Gehrmann-de Ridder, Thomas Gehrmann, Nigel Glover, Alexander Huss, Pier Francesco Monni, Emanuele Re, Luca Rottoli, Paolo Torrielli

► To cite this version:

Wojciech Bizoń, Xuan Chen, Aude Gehrmann-de Ridder, Thomas Gehrmann, Nigel Glover, et al. Fiducial distributions in Higgs and Drell-Yan production at N³LL+NNLO. Journal of High Energy Physics, 2018, 12, pp.132. 10.1007/JHEP12(2018)132. hal-01801919

HAL Id: hal-01801919

https://hal.science/hal-01801919v1

Submitted on 21 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Fiducial distributions in Higgs and Drell-Yan production at N³LL+NNLO

Wojciech Bizoń,ª Xuan Chen,ᵇ Aude Gehrmann-De Ridder,ᵇ,ᶜ Thomas Gehrmann,ᵇ

Nigel Glover,ᵈ Alexander Huss,ᵉ Pier Francesco Monni,ᵉ Emanuele Re,ᵉ,f

Luca Rottoliª and Paolo Torriellig

**ABSTRACT:** The perturbative description of certain differential distributions across a wide kinematic range requires the matching of fixed-order perturbation theory with resummation of large logarithmic corrections to all orders. We present precise matched predictions for transverse-momentum distributions in Higgs boson (H) and Drell-Yan pair (DY) production as well as for the closely related $\phi_\eta^*$ distribution at the LHC. The calculation is exclusive in the Born kinematics, and allows for arbitrary fiducial selection cuts on the decay products of the colour singlets, which is of primary relevance for experimental analyses. Our predictions feature very small residual scale uncertainties and display a good convergence of the perturbative series. A comparison of the predictions for DY observables to experimental data at 8 TeV shows a very good agreement within the quoted errors.

**KEYWORDS:** NLO Computations, QCD Phenomenology

**ARXIV ePRINT:** 1805.05916

<footnote>ªClarendon Laboratory, Rudolf Peierls Centre for Theoretical Physics, University of Oxford,
Parks Road, Oxford OX1 3PU, U.K.</footnote>

<footnote>ᵇDepartment of Physics, University of Zürich,
CH-8057 Zürich, Switzerland</footnote>

<footnote>ᶜInstitute for Theoretical Physics, ETH,
CH-8093 Zürich, Switzerland</footnote>

<footnote>ᵈInstitute for Particle Physics Phenomenology, Department of Physics, University of Durham,
Durham, DH1 3LE, U.K.</footnote>

<footnote>ᵉTheoretical Physics Department, CERN,
CH-1211 Geneva 23, Switzerland</footnote>

<footnote>f LAPTh, CNRS, Université Savoie Mont Blanc,
74940 Annecy, France</footnote>

<footnote n="9">⁹Dipartimento di Fisica and Arnold-Regge Center, Università di Torino,
and INFN, Sezione di Torino,
Via P. Giuria 1, I-10125 Torino, Italy
E-mail: wojciech.bizon@kit.edu, xuan.chen@uzh.ch, gehra@phys.ethz.ch,
thomas.gehrmann@uzh.ch, e.w.n.glover@durham.ac.uk,
alexander.huss@cern.ch, pier.monni@cern.ch,
emanuele.re@lapth.cnrs.fr, luca.rottoli@unimib.it,
torriell@to.infn.it</footnote>
# Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>2</td><td>Fixed order</td><td>4</td></tr><tr><td>3</td><td>Resummation</td><td>6</td></tr><tr><td>4</td><td>Matching to fixed order</td><td>10</td></tr><tr><td>4.1</td><td>Validation of the expansion and matching uncertainty for Higgs production</td><td>13</td></tr><tr><td>4.2</td><td>Validation of the expansion for Drell-Yan pair production</td><td>14</td></tr><tr><td>5</td><td>Results for Higgs production in HEFT</td><td>16</td></tr><tr><td>5.1</td><td>Matched predictions for inclusive Higgs</td><td>18</td></tr><tr><td>5.2</td><td>Matched predictions for fiducial H → γγ</td><td>19</td></tr><tr><td>6</td><td>Results for Drell-Yan production</td><td>21</td></tr><tr><td>6.1</td><td>Matched predictions for fiducial p<sub>t</sub><sup>2</sup> distributions</td><td>22</td></tr><tr><td>6.2</td><td>Matched predictions for fiducial φ<sub>η</sub><sup>*</sup> distributions</td><td>23</td></tr><tr><td>7</td><td>Conclusions</td><td>28</td></tr><tr><td>A</td><td>Formulae for the matching schemes</td><td>30</td></tr><tr><td>B</td><td>Formulae for N<sup>3</sup>LL resummation</td><td>32</td></tr></table>

## 1 Introduction

The accurate prediction and measurement of differential distributions is of primary importance for the LHC precision programme, especially in view of the absence of clear signals of new physics in the data collected so far. In this context, a special role is played by the kinematic distributions of a colour singlet produced in association with QCD radiation. These observables are often measured by reconstructing the decay products of the colour singlet (whenever possible), which are sensitive to the accompanying hadronic activity only through kinematic recoil. As a consequence, measurements of transverse and angular observables often lead to small experimental systematic uncertainties [1–8].

The implication of these precise measurements is twofold. On one hand, they can be used to fit the parameters of the SM Lagrangian (e.g. strong coupling constant, or masses) or to calibrate the models that typically enter the calculation of hadron-collider observables (like for instance collinear parton distribution functions (PDFs) [9], or non-perturbative corrections and transverse-momentum-dependent PDFs [10–12]). An example is given by
the differential distributions in Z- and W-boson production, that recently were exploited to perform very precise extractions of the W-boson mass [13] and to constrain the behaviour of some PDFs [14]. On the other hand, an excellent control over kinematic distributions is a way to set compelling constraints on new-physics models that would lead to mild shape distortions. An example is given by the sensitivity of the Higgs transverse-momentum ($p_t$) distribution to modification of the Yukawa couplings of the Higgs to quarks [15, 16].

In this article we present state-of-the-art predictions for a class of differential distributions both in Higgs boson (H) and Drell-Yan pair (DY) production. Specifically, we combine fixed-order calculations at next-to-next-to-leading order (NNLO) with the recently-obtained resummation of Sudakov logarithms to next-to-next-to-next-to-leading-logarithmic order (N³LL), for the transverse-momentum spectrum of the colour singlet, as well as for the angular variable $\phi_{\eta}^*$ [17]. In the following, for simplicity, we will collectively denote $p_t/M$ or $\phi_{\eta}^*$ by $v$, with $M$ representing the invariant mass of the colour singlet.

Inclusive and differential distributions for Higgs-boson production in gluon fusion are nowadays known with very high precision. The inclusive cross section has been computed to next-to-next-to-next-to-leading-order (N³LO) accuracy in QCD [18–24] in the heavy-top-quark limit. The impact of all-order effects due to a combined resummation of threshold and high-energy logarithms has been studied in detail, and at the current collider energies the corrections amount to a few-percent of the total cross section [25], indicating that the missing higher-order contributions are under good theoretical control. The state-of-the-art results for the Higgs transverse-momentum spectrum in fixed-order perturbation theory are the next-to-next-to-leading-order (NNLO) computations of refs. [26–29], which have been obtained in the heavy-top-quark limit. The effect of finite quark masses on differential distributions at next-to-leading order has been recently computed in refs. [30–35].

The state-of-the-art for the QCD corrections to differential distributions in DY production is at a similar level of accuracy. The total cross section is known fully differentially in the Born phase space up to NNLO [36–44], while differential distributions in transverse momentum were recently computed up to NNLO both for Z- [45–50] and W-boson [51–53] production. In the DY distributions, electroweak corrections become important especially at large transverse momenta, and they have been computed to NLO in [54–57].

Although fixed-order results are crucial to obtain reliable theoretical predictions away from the soft and collinear regions of the phase space ($v \sim 1$), it is well known that regions dominated by soft and collinear QCD radiation — which give rise to the bulk of the total cross section — are affected by large logarithmic terms of the form $\alpha_s^n \ln^k(1/v)/v$, with $k \le 2n - 1$, which spoil the convergence of the perturbative series at small $v$. In order to have a finite and well-behaved calculation in this limit, the subtraction of the infrared and collinear divergences requires an all-order resummation of the logarithmically divergent terms. The logarithmic accuracy is commonly defined in terms of the perturbative series of the logarithm of the cumulative cross section $\Sigma$ as

$$ \ln \Sigma(v) = \ln \int_0^v dv' \frac{d\Sigma(v')}{dv'} \\ = \sum_n \left\{ \mathcal{O}(\alpha_s^n \ln^{n+1}(1/v)) + \mathcal{O}(\alpha_s^n \ln^n(1/v)) + \mathcal{O}(\alpha_s^n \ln^{n-1}(1/v)) + \dots \right\}. \quad (1.1) $$
One refers to the dominant terms $\alpha_s^n \ln^{n+1}(1/v)$ as leading logarithmic (LL), to terms $\alpha_s^n \ln^n(1/v)$ as next-to-leading logarithmic (NLL), to $\alpha_s^n \ln^{n-1}(1/v)$ as next-to-next-to-leading logarithmic (NNLL), and so on.

The resummation of the $p_t$ spectrum of a heavy colour singlet is commonly performed in impact-parameter ($b$) space [58, 59], where the observable completely factorises and the resummed cross section takes an exponential form. Using the $b$-space formulation the Higgs $p_t$ spectrum was resummed at NNLL accuracy in refs. [60–62], following either the conventional approach of ref. [59], or a soft-collinear-effective-theory [63–66] (SCET) formulation of refs. [67, 68]. A study of the related theory uncertainties in the SCET formulation was presented in ref. [69]. In DY production, NNLL predictions for the transverse momentum of the color singlet as well as for $\phi_\eta^*$ were obtained in refs. [67, 70, 71]. The impact of both threshold and high-energy resummation on the small-transverse-momentum region was also studied in detail in refs. [72–80] and the effects were found to be quite moderate at LHC energies.

The problem of the resummation of the transverse-momentum distribution in direct ($p_t$) space received substantial attention throughout the years [81–83], but remained unsolved until recently. Due to the vectorial nature of $p_t$ (analogous considerations apply to $\phi_\eta^*$), it is indeed not possible to define a resummed cross section at a given logarithmic accuracy in direct space that is simultaneously free of both subleading-logarithmic contributions and spurious singularities at finite, non-zero values of $p_t$. A possible solution to the problem was recently proposed in refs. [84, 85], in whose formalism the resummation is performed by generating the relevant QCD radiation by means of a Monte Carlo (MC) algorithm. The resummation of the $p_t$ spectrum in momentum space has been also studied in ref. [86] within a SCET framework, where the renormalisation-group evolution is performed directly in $p_t$ space. An alternative technique to analytically transform the impact-parameter-space result into momentum space was recently proposed in ref. [87].

All the necessary ingredients for the $N^3LL$ resummation of $p_t$ (and $\phi_\eta^*$) spectra in color-singlet production have been computed in [88–93], and the four-loop cusp anomalous dimension has been recently obtained numerically in refs. [94, 95]. This has paved the way to more accurate theoretical results for transverse observables in the infrared region, like for instance the computation of the Higgs-transverse-momentum spectrum at $N^3LL$ matched to NNLO in refs. [85, 96]. In this manuscript, employing the direct-space resummation at $N^3LL$ accuracy of ref. [85] matched to NNLO, we present results for Higgs $p_t$ both at the inclusive level and with fiducial cuts on the decay products in the $H \to \gamma\gamma$ channel. We also consider Drell-Yan pair production and compute $N^3LL+NNLO$ predictions for the transverse momentum of the lepton pair and for the $\phi_\eta^*$ observable, comparing these results to ATLAS measurements at 8 TeV.

The article is organised as follows. In section 2 we discuss the computation of the NNLO differential distributions in DY and H production with the fixed-order parton-level code NNLOJET. Section 3 contains a brief review of the resummation for the $p_t$ and $\phi_\eta^*$ distributions using a momentum-space approach as implemented in the computer code RadISH, and in section 4 we discuss in detail the matching to fixed order together with the validation of our calculation. Section 5 reports the results for H production, while
the analogous results for DY production are reported in section 6. Section 7 contains our conclusions. We report the relevant formulae used for the matching in appendix A, while appendix B contains various quantities necessary for the resummation up to N³LL.

## 2 Fixed order

In this article we consider the production of either a Higgs boson or a leptonic Drell-Yan pair. In particular, the main focus lies in the description of the transverse-momentum spectrum and, in the case of DY production, of the closely related $\phi_{\eta}^*$ observable. These observables are studied in the context of matching the fixed-order calculation to a resummed prediction, and consequently the low- to intermediate-$p_t$ regimes are of particular interest.

For the Higgs production process, we therefore restrict ourselves to the region with $p_t^H \lesssim m_t$ where the HEFT description is appropriate. In this effective-field-theory framework, the top quark is integrated out in the large-top-mass limit ($m_t \to \infty$), giving rise to an effective operator that directly couples the Higgs field to the gluon field-strength tensor via [97–99]

$$ \mathcal{L}_{\text{HEFT}} = -\frac{\lambda}{4} G^{\mu\nu} G_{\mu\nu} H. \qquad (2.1) $$

The Wilson coefficient $\lambda$ is known to three-loop accuracy [100] and its renormalisation-scale dependence was studied in [29]. We consider the $p_t^H$ spectrum for both the inclusive production of an on-shell Higgs boson as well as including its decay into two photons. For the latter, the production and decay are treated in the narrow-width approximation and fiducial cuts, summarised in section 5, are applied on the photons in the final state.

For the DY process, we consider the full off-shell production of a charged lepton pair, including both the Z-boson and photon exchange contributions. Fiducial cuts are applied to the leptons in the final state and match the corresponding measurement performed by ATLAS at 8 TeV [101], which are summarised in section 6. We consider both the $p_t^Z$ spectrum as well as the $\phi_{\eta}^*$ distribution, which are further studied multi-differentially for different invariant-mass ($M_{\ell\ell}$) or rapidity ($Y_{\ell\ell}$) bins.

The differential distributions in $v = p_t/M$, $\phi_{\eta}^*$ for the production of a colour singlet at hadron colliders are indirectly generated through the recoil of the colour singlet against QCD radiation. The observables $v$ are therefore closely related to the $X + \text{jet}$ process with $X = H, Z$, where the jet requirement is replaced by a restriction on $v$ to be non-vanishing: $v \ge v_{\text{cut}} > 0$. The state-of-the-art fixed-order QCD predictions for this class of processes is at NNLO [26–29, 45–50]. Starting from the LO distributions, in which the colour singlet recoils against a single parton, the NNLO predictions receive contributions from configurations (with respect to LO) with two extra partons (RR: double-real corrections for H [102–104] and DY [105–109]), with one extra parton and one extra loop (RV: real-virtual corrections for H [110–112] and DY [105, 106, 113–116]) and with no extra parton but two extra loops (VV: double-virtual corrections for H [117] and DY [118–121]). Each of the three contributions is separately infrared divergent either in an implicit manner from phase-space regions where parton radiations become unresolved (soft and/or collinear) or in an explicit manner from divergent poles in virtual loop corrections. Only the sum of the three contributions is finite.
Our calculation is performed using the parton-level event generator NNLOJET, which implements the antenna subtraction method [122–124] to isolate infrared singularities and to enable their cancellation between different contributions prior to the numerical phase-space integration. The NNLO corrections for Higgs and DY production at finite $v$ are calculated using established implementations for $pp \to H + jet$ [29, 125] and $pp \to Z + jet$ [45–48] at NNLO, and it takes the schematic form:

$$ \sigma_{X+jet}^{\text{NNLO}} = \int_{\Phi_{X+3}} \left( d\sigma_{\text{NNLO}}^{\text{RR}} - d\sigma_{\text{NNLO}}^{S} \right) + \int_{\Phi_{X+2}} \left( d\sigma_{\text{NNLO}}^{\text{RV}} - d\sigma_{\text{NNLO}}^{T} \right) + \int_{\Phi_{X+1}} \left( d\sigma_{\text{NNLO}}^{\text{VV}} - d\sigma_{\text{NNLO}}^{U} \right). \quad (2.2) $$

The antenna subtraction terms, $d\sigma_{\text{NNLO}}^{S,T,U}$, for both Higgs and Drell-Yan related processes are constructed from antenna functions [126–131] to cancel infrared singularities between the contributions of different parton multiplicities. The integrals are performed over the phase space $\Phi_{X+1,2,3}$ corresponding to the production of the colour singlet in association with one, two or three partons in the final state. The integration of the final-state phase space is fully differential such that any infrared-safe observable $O$ can be studied through differential distributions as $d\sigma_{X+jet}^{\text{NNLO}}/dO$.

For large values of $v$ ($v \sim 1$), the phase-space integral in each line of eq. (2.2) is well defined and was calculated with high numerical precision in previous studies. Extending these predictions to smaller, but finite $v$ ($\sim 0.01$) becomes extremely challenging due to the wider dynamical range that is probed in the integration. Both the matrix elements and the subtraction terms grow rapidly in magnitude towards smaller values of $v$, thereby resulting in large numerical cancellations between them and rendering both the numerical precision and the stability of the results challenging. The finite remainder of such cancellations needs to be numerically stable in order to be consistently combined with resummed logarithmic corrections and extrapolated to the limit $v \to 0$. For this reason, the integration is performed separately for each individual initial-state partonic channel. We further split the integration region for each channel into multiple intervals in $v$, which are partially overlapping with each other. By carefully checking the consistency of the distributions in the overlapping region and using dedicated reweighting factors in each interval, we use NNLOJET to produce fixed-order predictions up to NNLO for values in $v$ down to $p_t = 2$ GeV and $\phi_{\eta}^* = 0.004$ [47].

The accuracy of the results obtained with the NNLOJET code for small $v$ has been systematically validated in ref. [96] by comparing fixed-order predictions of the Higgs boson transverse momentum distribution $d\sigma_{\text{NNLO}}/dp_t^H$ against the expansion of the $N^3\text{LL}$ resummation (obtained in the framework of soft-collinear effective field theory, SCET) to the respective order in the small $p_t^H$ region. This validation was performed for individual initial-state partonic channels down to $p_t^H = 0.7$ GeV.

As $v \to 0$, the final-state phase space $\Phi_{X+1,2,3}$ is reduced to the phase space of colour singlet production $\Phi_X$. The RR, RV, and VV contributions contain infrared di-
vergences with one extra unresolved parton that cannot be cancelled by the subtraction terms $d\sigma_{NNLO}^{S,T,U}$. These extra logarithmic divergences are cancelled by combining the fixed-order computation to a resummed calculation, where the logarithms in the fixed-order prediction are subtracted and replaced by a summation of the corresponding enhanced terms to all orders in perturbation theory. This operation is discussed in the next section, and more details on the combination of the two results are reported in appendix A.

## 3 Resummation

The approach developed in refs. [84, 85] uses the factorisation properties of the QCD squared amplitudes to devise a Monte Carlo formulation of the all-order calculation. In this framework, large logarithms are resummed directly in momentum space by effectively generating soft and/or collinear emissions in a fashion similar in spirit to an event generator.

To summarise the final result, we consider the cumulative distribution

$$ \Sigma(v) \equiv \int_0^v dv' \frac{d\Sigma(v')}{dv'} \qquad (3.1) $$

for an observable $v^{(l)} = V(\Phi_B, k_1, \dots, k_n)$, being either $p_t/M$ or $\phi_\eta^*$, in the presence of $n$ real emissions with momenta $k_1, \dots, k_n$. $\Sigma(v)$ can be expressed as

$$ \Sigma(v) = \int d\Phi_B \mathcal{V}(\Phi_B) \sum_{n=0}^{\infty} \int \prod_{i=1}^{n} [dk_i] |\mathcal{M}(\Phi_B, k_1, \dots, k_n)|^2 \Theta(v - V(\Phi_B, k_1, \dots, k_n)), \quad (3.2) $$

where $\mathcal{M}$ is the matrix element for $n$ real emissions and $\mathcal{V}(\Phi_B)$ denotes the resummed form factor that encodes the purely virtual corrections [132]. The phase spaces of the $i$-th emission $k_i$ and that of the Born configuration<fnref n="1" /> are denoted by $[dk_i]$ and $d\Phi_B$, respectively.

The recursive infrared and collinear (rIRC) safety [133] of the observable allows one to establish a well defined logarithmic counting in the squared amplitude [133, 134], and to systematically identify the contributions that enter at a given logarithmic order. In particular, the squared amplitude can be decomposed in terms of $n$-particle-correlated blocks, such that blocks with $n$ particles start contributing one logarithmic order higher than blocks with $n-1$ particles.

Eq. (3.2) contains exponentiated divergences of virtual origin in the $\mathcal{V}(\Phi_B)$ factor, as well as singularities in the real matrix elements, which appear at all perturbative orders. In order to handle such divergences, one can introduce a resolution scale $Q_0$ on the transverse momentum of the radiation: thanks to rIRC safety, unresolved real radiation (i.e. softer than $Q_0$) does not contribute to the observable's value, namely it can be neglected when computing $V(\Phi_B, k_1, \dots, k_n)$, thus it exponentiates and cancels the divergences contained in $\mathcal{V}(\Phi_B)$ at all orders. The precise definition of the unresolved radiation requires a careful clustering of momenta belonging to a given correlated block in order to be collinear safe. On

<footnote n="1">¹In the context of resummation, the Born configuration denotes the production of the colour-singlet state without any extra radiation. This should not be confused with the fixed-order counting of orders, where LO denotes the production of the colour-singlet state recoiling against a parton at finite transverse momentum.</footnote>
the other hand, the real radiation harder than the resolution scale (referred to as resolved) must be generated exclusively since it is constrained by the $\Theta$ function in eq. (3.2). rIRC safety also ensures that the dependence of the results upon $Q_0$ is power-like, hence the limit $Q_0 \to 0$ can be taken safely.

For observables which depend on the total transverse momentum of QCD radiation, such as $p_t$ or $\phi_\eta^*$, it is particularly convenient to set the resolution scale to a small fraction $\delta > 0$ of the transverse momentum of the block with largest $k_t$, hereby denoted by $\delta k_{t1}$, which allows for an efficient Monte Carlo implementation of the resulting resummed formula that can be used to simultaneously compute both $p_t$ and $\phi_\eta^*$.

Including terms up to N³LL, the cumulative cross section in momentum space can be recast in the following form [85]<fnref n="2" />

$$ \begin{aligned} \frac{d\Sigma(v)}{d\Phi_B} = & \int_0^\infty \frac{dk_{t1}}{k_{t1}} \mathcal{J}(k_{t1}) \frac{d\phi_1}{2\pi} \partial_{\tilde{L}} (-e^{-\tilde{R}(k_{t1})} \tilde{\mathcal{L}}^{\text{N}^3\text{LL}}(k_{t1})) \int d\mathcal{Z}[\{\tilde{R}', k_i\}] \Theta(v - V(\Phi_B, k_1, \dots, k_{n+1})) \\ & + \int_0^\infty \frac{dk_{t1}}{k_{t1}} \mathcal{J}(k_{t1}) \frac{d\phi_1}{2\pi} e^{-\tilde{R}(k_{t1})} \int d\mathcal{Z}[\{\tilde{R}', k_i\}] \int_0^1 \frac{d\zeta_s}{\zeta_s} \frac{d\phi_s}{2\pi} \left\{ \left( \tilde{R}'(k_{t1}) \tilde{\mathcal{L}}^{\text{NNLL}}(k_{t1}) - \partial_{\tilde{L}} \tilde{\mathcal{L}}^{\text{NNLL}}(k_{t1}) \right) \right. \\ & \left. \times \left( \tilde{R}''(k_{t1}) \ln \frac{1}{\zeta_s} + \frac{1}{2} \tilde{R}'''(k_{t1}) \ln^2 \frac{1}{\zeta_s} \right) - \tilde{R}'(k_{t1}) \left( \partial_{\tilde{L}} \tilde{\mathcal{L}}^{\text{NNLL}}(k_{t1}) - 2 \frac{\beta_0}{\pi} \alpha_s^2(k_{t1}) \hat{P}^{(0)} \otimes \tilde{\mathcal{L}}^{\text{NLL}}(k_{t1}) \ln \frac{1}{\zeta_s} \right) \right. \\ & \left. + \frac{\alpha_s^2(k_{t1})}{\pi^2} \hat{P}^{(0)} \otimes \hat{P}^{(0)} \otimes \tilde{\mathcal{L}}^{\text{NLL}}(k_{t1}) \right\} \left\{ \Theta(v - V(\Phi_B, k_1, \dots, k_{n+1}, k_s)) - \Theta(v - V(\Phi_B, k_1, \dots, k_{n+1})) \right\} \\ & + \frac{1}{2} \int_0^\infty \frac{dk_{t1}}{k_{t1}} \mathcal{J}(k_{t1}) \frac{d\phi_1}{2\pi} e^{-\tilde{R}(k_{t1})} \int d\mathcal{Z}[\{\tilde{R}', k_i\}] \int_0^1 \frac{d\zeta_{s1}}{\zeta_{s1}} \frac{d\phi_{s1}}{2\pi} \int_0^1 \frac{d\zeta_{s2}}{\zeta_{s2}} \frac{d\phi_{s2}}{2\pi} \tilde{R}'(k_{t1}) \\ & \times \left\{ \tilde{\mathcal{L}}^{\text{NLL}}(k_{t1}) (\tilde{R}''(k_{t1}))^2 \ln \frac{1}{\zeta_{s1}} \ln \frac{1}{\zeta_{s2}} - \partial_{\tilde{L}} \tilde{\mathcal{L}}^{\text{NLL}}(k_{t1}) \tilde{R}''(k_{t1}) \left( \ln \frac{1}{\zeta_{s1}} + \ln \frac{1}{\zeta_{s2}} \right) \right. \\ & \left. + \frac{\alpha_s^2(k_{t1})}{\pi^2} \hat{P}^{(0)} \otimes \hat{P}^{(0)} \otimes \tilde{\mathcal{L}}^{\text{NLL}}(k_{t1}) \right\} \\ & \times \left\{ \Theta(v - V(\Phi_B, k_1, \dots, k_{n+1}, k_{s1}, k_{s2})) - \Theta(v - V(\Phi_B, k_1, \dots, k_{n+1}, k_{s1})) - \Theta(v - V(\Phi_B, k_1, \dots, k_{n+1}, k_{s2})) + \Theta(v - V(\Phi_B, k_1, \dots, k_{n+1})) \right\} + \mathcal{O}\left(\alpha_s^n \ln^{2n-6} \frac{1}{v}\right), \end{aligned} \quad (3.3) $$

where $\zeta_{si} = k_{tsi}/k_{t1}$ and we introduced the notation $d\mathcal{Z}[\{\tilde{R}', k_i\}]$ to denote an ensemble that describes the emission of $n$ identical independent blocks [85]. The average of a function $G(\{\tilde{p}\}, \{k_i\})$ over the measure $d\mathcal{Z}$ is defined as $(\zeta_i = k_{ti}/k_{t1})$

$$ \int d\mathcal{Z}[\{\tilde{R}', k_i\}] G(\{\tilde{p}\}, \{k_i\}) = e^{-\tilde{R}'(k_{t1}) \ln \frac{1}{\delta}} \sum_{n=0}^{\infty} \frac{1}{n!} \prod_{i=2}^{n+1} \int_{\delta}^{1} \frac{d\zeta_i}{\zeta_i} \int_{0}^{2\pi} \frac{d\phi_i}{2\pi} \tilde{R}'(k_{t1}) G(\{\tilde{p}\}, k_1, \dots, k_{n+1}). \quad (3.4) $$

<footnote n="2">²We have split the result into a sum of three terms. The first term contains the full NLL corrections. The second term of eq. (3.3) (first set of curly brackets) starts contributing at NNLL accuracy, while the third term (second set of curly brackets) is purely N³LL.</footnote>
The $\ln 1/\delta$ divergence in the exponential prefactor of eq. (3.4) cancels exactly against that contained in the resolved real radiation, encoded in the nested sums of products on the right-hand side of the same equation. This ensures that the final result is therefore $\delta$-independent.

To obtain eq. (3.3) we used the fact that, for resolved radiation, $\zeta_i$ is a quantity of $O(1)$, which allows us to expand all ingredients in eq. (3.3) about $k_{t1}$, retaining only terms necessary for the desired logarithmic accuracy. We stress that this is allowed because of rIRC safety, which ensures that blocks with $k_{ti} \ll k_{t1}$ do not contribute to the value of the observable and are therefore fully cancelled by the term $\exp\{-\tilde{R}'(k_{t1}) \ln(1/\delta)\}$ of eq. (3.4). Although not strictly necessary, this expansion allows for a more efficient numerical implementation. The expansion gives rise to the terms $\tilde{R}^{(n)}$ which denote the derivatives of the radiator as

$$ \tilde{R}' = d\tilde{R}/d\tilde{L}, \quad \tilde{R}'' = d\tilde{R}'/d\tilde{L}, \quad \tilde{R}''' = d\tilde{R}''/d\tilde{L}, \qquad (3.5) $$

where $\tilde{R}$ takes the form

$$ \tilde{R}(k_{t1}) = -\tilde{L}g_1(\alpha_s\beta_0\tilde{L}) - g_2(\alpha_s\beta_0\tilde{L}) - \frac{\alpha_s}{\pi}g_3(\alpha_s\beta_0\tilde{L}) - \frac{\alpha_s^2}{\pi^2}g_4(\alpha_s\beta_0\tilde{L}), \qquad (3.6) $$

and $\alpha_s = \alpha_s(\mu_R)$. We report the functions $g_i$ in appendix B, and we refer to ref. [85] for further details. The function $g_4$ involves a contribution from the recently determined [95] four-loop cusp anomalous dimension $\Gamma_{\text{cusp}}^{(4)}$ that we report in eq. (B.12).

In previous $N^3\text{LL}$ resummation studies, $\Gamma_{\text{cusp}}^{(4)}$ was either neglected [85, 96] or extrapolated from its lower order contributions through a Padé approximation [135]. With the new result of [95] at hand, we could now explicitly verify that the numerical impact of $\Gamma_{\text{cusp}}^{(4)}$ is indeed very small (not visibly noticeable in the distributions), and well below other sources of parametric uncertainties that are discussed in the following.

The expression in eq. (3.3) would originally contain resummed logarithms of the form $\ln(Q/k_{t1})$, where $Q$ is the resummation scale, whose variation is used to probe the size of subleading logarithmic corrections not included in our result. In order to ensure that the resummation does not affect the hard region of the spectrum when matched to fixed order (see section 4), the resummed logarithms are supplemented with power-suppressed terms, negligible at small $k_{t1}$, that ensure resummation effects to vanish for $k_{t1} \gg Q$. Such modified logarithms $\tilde{L}$ are defined by constraining the rapidity integration of the real radiation to vanish at large transverse momenta. This is done by mapping the limit $k_{t1} \to Q$ onto $k_{t1} \to \infty$ in all terms of eq. (3.3), with the exception of the observable's measurement function. A convenient choice of such a mapping is

$$ \ln \frac{Q}{k_{t1}} \to \tilde{L} = \frac{1}{p} \ln \left( \left( \frac{Q}{k_{t1}} \right)^p + 1 \right), \qquad (3.7) $$

where $p$ is a positive real parameter chosen in such a way that the resummed differential distribution vanishes faster than the fixed-order one at large $v$, with slope $(1/v)^{p+1}$. The above prescription comes with the prefactor $\mathcal{J}$, defined as

$$ \mathcal{J}(k_{t1}) = \left( \frac{Q}{k_{t1}} \right)^p \left( 1 + \left( \frac{Q}{k_{t1}} \right)^p \right)^{-1}. \qquad (3.8) $$
This corresponds to the Jacobian for the transformation (3.7), and ensures the absence of fractional (although power suppressed) $\alpha_s$ powers in the final distribution [85]. This factor, once again, leaves the small $k_{t1}$ region untouched, and only modifies the large $p_t$ region by power-suppressed effects. Although this procedure seems a simple change of variables, we stress that the observable's measurement function (i.e. the $\Theta$ function in eq. (3.3)) is not affected by this prescription. As a consequence, the final result will depend on the parameter $p$ through power-suppressed terms.

The factors $\tilde{\mathcal{L}}$ contain the parton luminosities up to N³LL, multiplied by the Born-level squared, and virtual amplitudes. They are defined as

$$ \tilde{\mathcal{L}}_{\text{NLL}}(k_{t1}) = \sum_{c,c'} \frac{d|M_B|_{cc'}^2}{d\Phi_B} f_c(\mu_F e^{-\tilde{L}}, x_1) f_{c'}(\mu_F e^{-\tilde{L}}, x_2), \quad (3.9) $$

$$ \tilde{\mathcal{L}}_{\text{NNLL}}(k_{t1}) = \sum_{c,c'} \frac{d|M_B|_{cc'}^2}{d\Phi_B} \sum_{i,j} \int_{x_1}^{1} \frac{dz_1}{z_1} \int_{x_2}^{1} \frac{dz_2}{z_2} f_i(\mu_F e^{-\tilde{L}}, \frac{x_1}{z_1}) f_j(\mu_F e^{-\tilde{L}}, \frac{x_2}{z_2}) \times \left\{ \delta_{ci}\delta_{c'j}\delta(1-z_1)\delta(1-z_2) \left(1 + \frac{\alpha_s(\mu_R)}{2\pi} \tilde{H}^{(1)}(\mu_R, x_Q)\right) + \frac{\alpha_s(\mu_R)}{2\pi} \frac{1}{1-2\alpha_s(\mu_R)\beta_0 \tilde{L}} \left(\tilde{C}_{ci}^{(1)}(z_1, \mu_F, x_Q)\delta(1-z_2)\delta_{c'j} + \{z_1 \leftrightarrow z_2; c, i \leftrightarrow c', j\}\right) \right\}, \quad (3.10) $$

$$ \tilde{\mathcal{L}}_{\text{N}^3\text{LL}}(k_{t1}) = \sum_{c,c'} \frac{d|M_B|_{cc'}^2}{d\Phi_B} \sum_{i,j} \int_{x_1}^{1} \frac{dz_1}{z_1} \int_{x_2}^{1} \frac{dz_2}{z_2} f_i(\mu_F e^{-\tilde{L}}, \frac{x_1}{z_1}) f_j(\mu_F e^{-\tilde{L}}, \frac{x_2}{z_2}) \times \left\{ \begin{aligned} & \delta_{ci}\delta_{c'j}\delta(1-z_1)\delta(1-z_2) \left(1 + \frac{\alpha_s(\mu_R)}{2\pi} \tilde{H}^{(1)}(\mu_R, x_Q) + \frac{\alpha_s^2(\mu_R)}{(2\pi)^2} \tilde{H}^{(2)}(\mu_R, x_Q)\right) \\ & + \frac{\alpha_s(\mu_R)}{2\pi} \frac{1}{1-2\alpha_s(\mu_R)\beta_0 \tilde{L}} \left(1 - \alpha_s(\mu_R) \frac{\beta_1}{\beta_0} \frac{\ln(1-2\alpha_s(\mu_R)\beta_0 \tilde{L})}{1-2\alpha_s(\mu_R)\beta_0 \tilde{L}}\right) \\ & \times \left(\tilde{C}_{ci}^{(1)}(z_1, \mu_F, x_Q)\delta(1-z_2)\delta_{c'j} + \{z_1 \leftrightarrow z_2; c, i \leftrightarrow c', j\}\right) \\ & + \frac{\alpha_s^2(\mu_R)}{(2\pi)^2} \frac{1}{(1-2\alpha_s(\mu_R)\beta_0 \tilde{L})^2} \left(\tilde{C}_{ci}^{(2)}(z_1, \mu_F, x_Q)\delta(1-z_2)\delta_{c'j} + \{z_1 \leftrightarrow z_2; c, i \leftrightarrow c', j\}\right) \\ & + \frac{\alpha_s^2(\mu_R)}{(2\pi)^2} \frac{1}{(1-2\alpha_s(\mu_R)\beta_0 \tilde{L})^2} \left(\tilde{C}_{ci}^{(1)}(z_1, \mu_F, x_Q)\tilde{C}_{c'j}^{(1)}(z_2, \mu_F, x_Q) + G_{ci}^{(1)}(z_1)G_{c'j}^{(1)}(z_2)\right) \\ & + \frac{\alpha_s^2(\mu_R)}{(2\pi)^2} \tilde{H}^{(1)}(\mu_R, x_Q) \frac{1}{1-2\alpha_s(\mu_R)\beta_0 \tilde{L}} \left(\tilde{C}_{ci}^{(1)}(z_1, \mu_F, x_Q)\delta(1-z_2)\delta_{c'j} + \{z_1 \leftrightarrow z_2; c, i \leftrightarrow c', j\}\right) \end{aligned} \right\}. \quad (3.11) $$

where

$$ x_1 = \frac{M}{\sqrt{s}} e^Y, \qquad x_2 = \frac{M}{\sqrt{s}} e^{-Y}, \quad (3.12) $$

$Y$ is the rapidity of the colour singlet in the centre-of-mass frame of the collision at the Born-level, $|M_B|_{cc'}^2$ is the Born-level squared matrix element, and $x_Q = Q/M$. The above luminosities contain the NLO and NNLO coefficient functions $\tilde{C}_{ci}^{(n)}$ for Higgs and Drell-Yan production [88–91], as well as the hard virtual corrections $\tilde{H}^{(n)}$. A precise definition is given in section 4 of ref. [85], and the relevant formulae are also reported in appendix B.
Finally, we define the convolution of a regularised splitting function $\hat{P}$ [136, 137] with the coefficient $\tilde{\mathcal{L}}_{\text{NLL}}$ as

$$ \hat{P}^{(0)} \otimes \tilde{\mathcal{L}}_{\text{NLL}}(k_{t1}) = \sum_{c,c'} \frac{d|M_B|_{cc'}^2}{d\Phi_B} \left\{ \begin{aligned} & (\hat{P}^{(0)} \otimes f)_c (\mu_F e^{-\tilde{L}}, x_1) f_{c'} (\mu_F e^{-\tilde{L}}, x_2) \\ & + f_c (\mu_F e^{-\tilde{L}}, x_1) (\hat{P}^{(0)} \otimes f)_{c'} (\mu_F e^{-\tilde{L}}, x_2) \end{aligned} \right\}. \quad (3.13) $$

The term $\hat{P}^{(0)} \otimes \hat{P}^{(0)} \otimes \tilde{\mathcal{L}}_{\text{NLL}}(k_{t1})$ is to be interpreted in the same way. Moreover, the explicit factors of the strong coupling evaluated at $k_{t1}$ in eq. (3.3) are defined as

$$ \alpha_s(k_{t1}) = \frac{\alpha_s(\mu_R)}{1 - 2\alpha_s(\mu_R)\beta_0 \tilde{L}}. \quad (3.14) $$

## 4 Matching to fixed order

In this section we discuss the matching of the resummed and the fixed-order results. Since we work at the level of the cumulative distribution $\Sigma$, we define the analogue of eq. (3.1) for the fixed-order prediction as

$$ \Sigma^{\text{N}^3\text{LO}}(v) = \sigma_{\text{tot}}^{\text{N}^3\text{LO}} - \int_v^\infty dv' \frac{d\Sigma^{\text{NNLO}}(v')}{dv'}, \quad (4.1) $$

where $\sigma_{\text{tot}}^{\text{N}^3\text{LO}}$ is the total cross section for the considered processes and $d\Sigma^{\text{NNLO}}/dv'$ denotes the NNLO differential distribution.

For inclusive Higgs production, the transverse-momentum distribution at NNLO was obtained in refs. [26–29], while the N$^3$LO total cross section has been computed in refs. [23, 24]. On the other hand, the N$^3$LO cross section within fiducial cuts on the Born kinematics is currently unknown. Since in this article we address differential distributions for $H \to \gamma\gamma$ with fiducial cuts, we approximate the N$^3$LO correction to $\sigma_{\text{tot}}^{\text{N}^3\text{LO}}$ by rescaling the NNLO fiducial cross section by the inclusive (i.e. without fiducial cuts) N$^3$LO/NNLO K factor. We stress that, at the level of the differential distributions we are interested in, this approximation is formally a N$^4$LL effect, and it lies beyond the accuracy considered in this study.

For DY production, the differential distributions to NNLO were obtained in refs. [47, 49]. We set to zero the unknown N$^3$LO correction to the total cross section, observing once again that its contribution to the distributions derived here is subleading.

In order to assess the uncertainty associated with the matching procedure, we consider here two different matching schemes. The first scheme we introduce is the common additive scheme defined as

$$ \Sigma_{\text{add}}^{\text{MAT}}(v) = \Sigma^{\text{N}^3\text{LL}}(v) + \Sigma^{\text{N}^3\text{LO}}(v) - \Sigma^{\text{EXP}}(v), \quad (4.2) $$

where $\Sigma^{\text{EXP}}$ denotes the expansion of the resummation formula $\Sigma^{\text{N}^3\text{LL}}$ to N$^3$LO.
The second scheme we consider belongs to the class of multiplicative schemes similar to those defined in refs. [138–140], and it is schematically defined as

$$ \Sigma_{\text{mult}}^{\text{MAT}}(v) = \Sigma^{\text{N}^3\text{LL}}(v) \left[ \frac{\Sigma^{\text{N}^3\text{LO}}(v)}{\Sigma^{\text{EXP}}(v)} \right]_{\text{EXPANDED TO N}^3\text{LO}}, \quad (4.3) $$

where the quantity in square brackets is expanded to N³LO. The two schemes (4.2), (4.3) are equivalent at the perturbative order we are working at, and differ by N⁴LO and N⁴LL terms. The main difference between the two schemes is that in the multiplicative approach, unlike in the additive one, higher-order corrections are damped by the resummation factor ΣN³LL at low v. Moreover, this damping occurs in the region where the fixed-order result may be occasionally affected by numerical instabilities, hence allowing for a stable matched distribution even with limited statistics for the NNLO component.

One advantage of the multiplicative solution is that the N³LO constant terms, of formal N⁴LL accuracy, are automatically extracted from the fixed order in the matching procedure, whenever the N³LO total cross section is known. We recall that eq. (3.3) resums all towers of ln(1/v) up to N³LL, defined at the level of the logarithm of Σ (1.1). At this order, one predicts correctly all logarithmic terms up to, and including, α_s^n ln(2n-5)(1/v) in the expanded formula for Σ, while terms of order α_s^n ln(2n-6)(1/v) would be modified by including N⁴LL corrections.

The inclusion of constant terms of order O(α_s<fnref n="3" />) relative to Born level in the resummed formula, of formal N⁴LL accuracy, extends the prediction to all terms of order α_s^n ln(2n-6)(1/v) in the expanded formula for Σ. Indeed these terms, which contain the N³LO collinear coefficient functions and three-loop virtual corrections, would multiply the Sudakov e⁻R(k_t1) in the resummed formula (3.3) starting at N⁴LL. Since they are currently unknown analytically, in an additive matching these terms are simply added to the resummed cumulative result, and disappear at the level of the differential distribution. On the other hand, in a multiplicative scheme, they multiply the resummed cross section and hence correctly include a whole new tower of N⁴LL terms α_s^n ln(2n-6)(1/v) in the expanded formula for the matched cumulative cross section ΣMAT.<fnref n="3" /> We stress that this, as pointed out above, requires the knowledge of the N³LO cross section in the considered fiducial volume. This is currently only known in the case of fully inclusive Higgs production, whose results are presented in section 5.1. In the remaining studies of fiducial distributions, both for Higgs in section 5.2, and for DY in section 6, the N³LO cross sections are approximated, as described at the beginning of this section, and hence the tower of N⁴LL terms α_s^n ln(2n-6)(1/v) in Σ is not fully included.

However, there is a drawback in using eq. (4.3) as is. Indeed, in the limit $\tilde{L} \to 0$, ΣN³LL tends to the integral of $\tilde{L}_{N^3LL}(\mu_F)$ (defined in eq. (3.11)) over $\Phi_B$, evaluated at $\tilde{L} = 0$. Therefore, the fixed-order result ΣN³LO at large v receives a spurious correction of relative order α_s⁴

$$ \Sigma_{\text{mult}}^{\text{MAT}}(v) \sim \Sigma^{\text{N}^3\text{LO}}(v) (1 + \mathcal{O}(\alpha_s^4)). \quad (4.4) $$

<footnote n="3">³Notice that this does not imply that the whole class of N⁴LL terms is included. This would instead require all terms of the form α_s^n ln(n⁻³)(1/v) in lnΣ, eq. (1.1), which would predict correctly all terms α_s^n ln(2n-6)(1/v) and α_s^n ln(2n-7)(1/v) in the expanded Σ.</footnote>
Despite being formally of higher order, this effect can be moderately sizeable in processes with large K factors, such as Higgs production. There are different possible solutions to this problem. In ref. [85] the resummed component (as well as the relative expansion) was modified by introducing a damping factor as

$$ \Sigma^{\mathrm{N}^3\mathrm{LL}} \rightarrow (\Sigma^{\mathrm{N}^3\mathrm{LL}})^Z, \qquad (4.5) $$

where Z is a v-dependent exponent that effectively acts as a smoothened Θ function that tends to zero at large v. This solution, however, introduces new parameters that control the scaling of the damping factor Z (see section 4.2 of ref. [85] for details). In this article we adopt a simpler solution, which avoids the introduction of extra parameters in the matching scheme. To this end, we define the multiplicative matching scheme by normalising the resummed prefactor to its asymptotic $\tilde{L} \to 0$ value. This is simply given by the integral over the Born phase space $\Phi_B$ of the $\tilde{L} \to 0$ limit of $\tilde{\mathcal{L}}_{\mathrm{N}^3\mathrm{LL}}$ (that we report in eq. (A.5))

$$ \Sigma_{\mathrm{asym.}}^{\mathrm{N}^3\mathrm{LL}} = \int_{\mathrm{with cuts}} d\Phi_B \left( \lim_{\tilde{L} \to 0} \tilde{\mathcal{L}}_{\mathrm{N}^3\mathrm{LL}} \right), \qquad (4.6) $$

where the integration over $\Phi_B$ is performed by taking into account the phase-space cuts of the experimental analysis.

We thus obtain

$$ \Sigma_{\mathrm{mult}}^{\mathrm{MAT}}(v) = \frac{\Sigma^{\mathrm{N}^3\mathrm{LL}}(v)}{\Sigma_{\mathrm{asym.}}^{\mathrm{N}^3\mathrm{LL}}} \left[ \Sigma_{\mathrm{asym.}}^{\mathrm{N}^3\mathrm{LL}} \frac{\Sigma^{\mathrm{N}^3\mathrm{LO}}(v)}{\Sigma^{\mathrm{EXP}}(v)} \right]_{\mathrm{EXPANDED TO N}^3\mathrm{LO}}, \qquad (4.7) $$

where

$$ \Sigma^{\mathrm{N}^3\mathrm{LL}}(v) \xrightarrow{v \gg Q/M} \Sigma_{\mathrm{asym.}}^{\mathrm{N}^3\mathrm{LL}} \qquad (4.8) $$

and the whole squared bracket in eq. (4.7) is expanded to $\mathrm{N}^3\mathrm{LO}$. This ensures that, in the $v \gg Q/M$ limit, eq. (4.7) reproduces by construction the fixed-order result, and no large spurious higher-order corrections arise in this region. The detailed matching formulae for the two schemes considered in our analysis are reported in appendix A.

In order to estimate the systematic uncertainty associated with the choice of the matching scheme, a consistent comparison between the two will be performed in the next section considering inclusive Higgs production as a case study.

Before we proceed with the results, we stress that in the remainder of this article we will only focus on differential distributions rather than on cumulative ones. Therefore, at the level of the spectrum, in our notation we will drop one order in the fixed-order counting, so that the derivative of $\Sigma^{\mathrm{N}^3\mathrm{LO}}$ will be referred to as a NNLO distribution, and analogously for the lower-order cases.

In the next two subsections we perform some validation studies both for Higgs (section 4.1) and DY (section 4.2) production, where we compare the fixed-order calculation in the deep infrared regime to the expansion of the resummed result. Moreover, we discuss the uncertainty associated with the choice of the matching scheme, and estimate it through a comparison of the two prescriptions defined above for a case study.
## 4.1 Validation of the expansion and matching uncertainty for Higgs production

To perform the matching to fixed order, the resummation formula (3.3) is expanded up to the third order in the strong coupling. To obtain the expanded results, one can directly set the resolution scale $\delta$ to zero, since the cancellation of IRC divergences is manifest. In figure 1 we show the comparison between the expansion of the N³LL resummed cross section and the fixed order for the differential distribution of $p_t^H$ both at NLO (left plot) and at NNLO (right plot). We remind the reader that at the level of the differential distribution NNLO denotes the derivative of the N³LO cumulant, and similarly for lower orders.

In figure 1 we see that below $p_t^H \sim 10$ GeV the fixed-order and the expansion of the resummation are in excellent agreement, and that the size of non-logarithmic terms in the perturbative series remains moderate up to $p_t^H \sim 50$ GeV.

It is instructive to further investigate the difference between the fixed order and the expansion of the resummation formula in the region of very small $p_t^H$. In particular, we consider the differential distribution

$$ \frac{d\Sigma(p_t^H)}{d\ln(p_t^H/\text{GeV})}, \qquad (4.9) $$

in order to highlight potential logarithmic differences in the $p_t^H \to 0$ region. A similar validation of the NNLO $p_t^H$ distribution has been performed in ref. [96]. The result of our comparison is displayed in the left panel of figure 2. The dashed green line shows the difference between the NNLO distribution and the $O(\alpha_s^3)$ expansion of the NNLL resummation. As one expects, at small $p_t^H$ the two predictions for the cumulative distribution differ by a double-logarithmic term (due to the absence of the NNLO coefficient functions and of the two-loop virtual corrections in the NNLL result), which induces a linear slope at the level of the differential distribution (4.9). When we include the N³LL corrections (solid red line), the difference between the two curves tends to zero, hence proving the consistency between the two predictions. For comparison, the difference between the NLO and NNLL (cyan dot-dashed line) is also reported. The right panel of figure 2 shows the difference between the NNLO coefficient and the corresponding expansion of the N³LL resummation at the same order. The lower inset of the same figure shows the ratio of the above difference to the NNLO coefficient, which helps quantify the relative difference.

As a check on the theoretical setup that will be used in the next sections, it is interesting to compare the predictions for the $p_t^H$ spectrum obtained with the two matching schemes defined in eqs. (4.2) and (4.7). In order to compare the multiplicative and additive schemes on an equal footing, hence including the same ingredients for both schemes, in this section we consider a matching to NNLO at the level of the cumulative cross section that will allow us to estimate the systematic uncertainty associated with the choice of the matching scheme. In this case the resummed cross section is defined as in eqs. (4.2) and (4.7) with the obvious replacement of N³LO by NNLO. The result of the comparison is reported in figure 3. We observe a very good agreement between the two matching schemes, which is a sign of robustness of the predictions shown below. The lower panel of figure 3 shows the relative uncertainty bands obtained within the two schemes, where each prediction is
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Comparison between the fixed-order transverse momentum distribution for Higgs boson production at $\sqrt{s} = 13$ TeV at NLO (left) and NNLO (right) and the expansion of the N³LL resummation formula given in eq. (3.3) to the corresponding order, i.e. $O(\alpha_s^4)$ and $O(\alpha_s^5)$ (namely $O(\alpha_s^2)$ and $O(\alpha_s^3)$ relative to Born), respectively.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Left: difference between the full NLO and NNLO $p_t^H$ distribution and the expansion of the NNLL and N³LL resummation formulae (3.3) to the respective perturbative order. Right: difference between the fixed-order NNLO coefficient, i.e. the $O(\alpha_s^5)$ term alone, and the corresponding coefficient obtained from the expansion of the N³LL resummation.</figcaption></figure>

divided by its own central value. The theory uncertainties have a very similar pattern. Given that the difference between the two schemes is always quite moderate with respect to the scale uncertainty, in the following we decide to proceed with the multiplicative prescription (4.7) as our default. We find analogous conclusions for DY production, and therefore we choose not to report this further comparison here.

## 4.2 Validation of the expansion for Drell-Yan pair production

Similarly to the validation performed for inclusive Higgs production, in this section we consider the difference between the NNLO differential distribution and the corresponding expansion of the N³LL resummed calculation. In particular, we focus on the differential
<figure><img src="image_4.png" /><figcaption>Figure 3. Comparison between additive and multiplicative matching schemes at N³LL+NLO for the transverse momentum distribution for Higgs boson production at $\sqrt{s}$ = 13 TeV. The lower panel shows the relative uncertainty bands obtained within the two schemes.</figcaption></figure>

distribution

$$ \frac{d\Sigma(p_t^Z)}{d\ln(p_t^Z/\text{GeV})}, \qquad (4.10) $$

in order to highlight potential logarithmic differences in the $p_t^Z \to 0$ region.

To perform the validation we consider 8 TeV pp collisions with NNPDF3.0 parton densities [141], and we work within an inclusive setup requiring

$$ 80 \text{ GeV} < M_{\ell\ell} < 100 \text{ GeV}, \qquad (4.11) $$

and setting the scales to $\mu_R = \mu_F = M_Z$ with $x_Q = Q/M_{\ell\ell} = 1$. This inclusive setup is chosen as to avoid any potential complications due to the use of fiducial cuts, as well as dynamical scales, that act differently on the fixed-order and resummed calculations. Indeed, at variance with the case of the fixed-order calculation, in the resummation both fiducial cuts and dynamical scales are always defined at level of the Born (i.e. $Z + 0$ jet) phase space, which differs from the definition used in the fixed-order calculation unless the extra QCD radiation is extremely soft or collinear to the beam. As a consequence, employing fiducial cuts and/or dynamical scales may necessitate going to smaller values of $p_t^Z$ in order to see a convergence of the fixed-order to the expansion of the resummation.

The results of the comparison are shown in figure 4. The left panel displays the difference between the NLO distribution and the expansion of the NNLL resummation to second order (cyan dot-dashed line), and between the NNLO distribution and the expansion of the N³LL resummation to third order (solid red line). In both cases one expects the differences to approach zero at small $p_t^Z$, which is well confirmed by the plot. In addition, we report on the difference between the NNLO distribution and the expansion of the NNLL resummation to third order given by the dashed green line. Due to missing double-logarithmic terms in the NNLL expansion, a non-vanishing slope is expected in the low-$p_t^Z$ region, which is suggested by the green curve within statical uncertainties. In order to single out the contribution of the NNLO correction, in the right panel of figure 4 we show
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Left: difference between the full NLO and NNLO $p_t^Z$ distribution and the expansion of the NNLL and N$^3$LL resummation formulae (3.3) to the respective perturbative order. Right: difference between the fixed-order NNLO coefficient, i.e. the $O(\alpha_s^3)$ term alone, and the corresponding coefficient obtained from the expansion of the N$^3$LL resummation.</figcaption></figure>

the difference between the NNLO coefficient alone, and the corresponding coefficient in the expansion of the N$^3$LL resummation. As expected, such a difference asymptotically tends to zero for small $p_t^Z$ values.

In addition to the validation of the full $p_t^Z$ spectrum shown in figure 4, we have further performed the analogous checks for the individual partonic channels which are summarised in figure 5. To this end, we have computed the fixed-order NNLO contribution to the $p_t^Z$ distribution down to $p_t^Z \sim 0.5$ GeV with uncertainties at the 10% level. We can clearly observe that the fixed-order prediction is in excellent agreement with the prediction from the resummed calculation for all partonic configurations. The respective bottom panels in each figure show the difference between the two predictions, which for all channels approach zero in the limit $p_t^Z \to 0$. This is an excellent cross-check of the two calculations, which proves the good numerical stability of the NNLO distributions down to the deep infrared regime.

# 5 Results for Higgs production in HEFT

In this section we present our predictions for the $p_t^H$ spectrum both in inclusive $pp \to H$ production, and in the $pp \to H \to \gamma\gamma$ channel with fiducial cuts. The computational setup is the same for both analyses, and all results presented below are obtained in the heavy-top-quark limit. We consider collisions at 13 TeV, and use parton densities from the PDF4LHC15_nnlo_mc set [141–146]. The value of the parameter $p$ appearing in the definition of the modified logarithms $\tilde{L}$ is chosen considering the scaling of the spectrum in the hard region, so as to make the matching to the fixed order smooth there. We set $p = 4$ as our reference value, but nevertheless have checked that a variation of $p$ by one unit does not induce any significant differences.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Validation between the fixed-oder coefficients (at NLO and NNLO) and the corresponding expansion of the resummed prediction (at NNLL and N³LL) for the individual partonic channels, with $L = \ln(p_t^Z/\text{GeV})$. Note that in contrast to figure 4, the curves labelled as "NNLL" only comprises term of $O(\alpha_s^2)$ and does not include higher-order $O(\alpha_s^3)$ terms.</figcaption></figure>

We set the central renormalisation and factorisation scales as $\mu_R = \mu_F = m_H/2$, with $m_H = 125$ GeV, while the resummation scale is chosen to be $x_Q = Q/m_H = 1/2$. We estimate the perturbative uncertainty by performing a seven-scale variation of $\mu_R$, $\mu_F$ by a factor of two in either direction, while keeping $1/2 < \mu_R/\mu_F < 2$ and $x_Q = 1/2$; Moreover,
for central $\mu_R$ and $\mu_F$ scales, $x_Q$ is varied around its central value by a factor of two. The quoted theoretical error is defined as the envelope of all the above variations. We discuss the results for inclusive production in section 5.1, and then present the predictions for the fiducial distributions in section 5.2.

## 5.1 Matched predictions for inclusive Higgs

We start by quantifying the size of the N$^3$LL effects compared to NNLL resummation. In the left plot of figure 6 we compare the differential distributions at N$^3$LL+NLO and NNLL+NLO in the small-$p_t^H$ region. The lower panel of the plot shows the ratio of both predictions to the central line of the N$^3$LL+NLO band, which corresponds to central scales in our setup. We observe that N$^3$LL corrections are very moderate in size, with effects of order 2% on the central prediction in most of the displayed range, growing up to at most 5% only in the region of extremely low $p_t^H$. The central N$^3$LL+NLO result is entirely contained in the NNLL+NLO uncertainty band. On the other hand, the inclusion of the N$^3$LL corrections reduces the perturbative uncertainty for $p_t^H \lesssim 5$ GeV.

The right plot of figure 6 shows the same comparison for the matching to NNLO. The effect of the N$^3$LL corrections is consistent with the previous order, with a percent-level correction in most of the range, growing up to 5% at very small $p_t^H$. Similarly, the perturbative uncertainty is significantly reduced below 10 GeV with respect to the NNLL+NNLO case. It is important to stress that in the NNLL+NNLO matching the fixed order and the expansion of the resummation differ by a divergent term $\sim 1/p_t^H$ at small $p_t^H$. The fact that the divergence is not visible in the distribution reported in the upper panel of figure 6 is entirely due to the nature of the multiplicative scheme, which ensures that the distribution follows the resummation scaling at small $p_t^H$, therefore damping the divergence. A multiplicative matching of N$^3$LL resummation to NNLO was already shown in ref. [85], where however no significant reduction in the uncertainty band at small $p_t^H$ was observed in that case. This feature was due to the limited statistics of the fixed-order distributions used in that analysis at small $p_t^H$, whose fluctuations dominated the uncertainty band at very small transverse momentum. An additive matching of N$^3$LL to NNLO was recently performed in ref. [96].

Next, we consider the comparison between the matched prediction and the fixed-order one. Figure 7 shows this comparison for two different central scales. The left plot is obtained with central $\mu_F = \mu_R = m_H/2$, while the right plot is obtained with $\mu_F = \mu_R = m_H$. The rest of the setup is kept as described above. We observe that at $\mu_F = \mu_R = m_H/2$ the uncertainty band is affected by cancellations in the scale variation, which accidentally lead to a small perturbative uncertainty. Choosing $m_H$ as a central scale (right plot of figure 7) leads to a broader uncertainty band resulting in a more robust estimate of the perturbative error. This is particularly the case for predictions above 50 GeV, where resummation effects are progressively less important. We notice indeed that in both cases the effect of resummation starts to be increasingly relevant for $p_t^H \lesssim 40$ GeV.

In the following we choose $m_H/2$ as a central scale. Nevertheless, we stress that a comparison to data (not performed here for Higgs boson production) will require a study of different central-scale choices.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Comparison between different combinations of fixed-order (NLO and NNLO) and resummation (NNLL and N³LL) for the transverse momentum distribution for Higgs boson production at √s = 13 TeV. Left: NLO and Right: NNLO. The lower panel shows the ratio of predictions to that obtained with N³LL resummation.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Comparison of the transverse momentum distribution for Higgs boson production at NNLO and N³LL+NNLO for a central scale choice of μR = μF = mH/2 (left) and μR = μF = mH (right). In both cases, Q = mH/2. The lower panel shows the ratio to the N³LL+NNLO prediction.</figcaption></figure>

To conclude, figure 8 reports the comparison between our best prediction (N³LL+NNLO), the NNLL+NLO, and the NNLO distributions. The plot shows a very good convergence of the predictions at different perturbative orders, with a significant reduction of the scale uncertainty in the whole kinematic range considered here.

## 5.2 Matched predictions for fiducial $H \to \gamma\gamma$

Experimental measurements are performed within a fiducial phase-space volume, defined in order to comply with the detector geometry and to enhance signal sensitivity. On the theoretical side it is therefore highly desirable to provide predictions that exactly match the experimental setup. The availability of matched predictions that are fully differential in
<figure><img src="image_9.png" /><figcaption>Figure 8. Comparison of the transverse momentum distribution for Higgs boson production between N³LL+NNLO, NNLL+NLO, and NNLO at central scale choice of $\mu_R = \mu_F = m_H/2$. The lower panel shows the ratio to the N³LL+NNLO prediction.</figcaption></figure>

the Born phase space also allows for a direct comparison to data without relying on Monte Carlo modeling of acceptances. In this section we consider the process $pp \to H \to \gamma\gamma$ and, in particular, we focus on the transverse momentum of the $\gamma\gamma$ system in the presence of fiducial cuts.

The fiducial volume is defined by the set of cuts detailed below [7]

$$ \begin{gathered} \min(p_t^{\gamma 1}, p_t^{\gamma 2}) > 31.25 \text{ GeV}, \qquad \max(p_t^{\gamma 1}, p_t^{\gamma 2}) > 43.75 \text{ GeV}, \\ 0 < |\eta^{\gamma 1,2}| < 1.37 \text{ or } 1.52 < |\eta^{\gamma 1,2}| < 2.37, \qquad |Y_{\gamma\gamma}| < 2.37, \end{gathered} \quad (5.1) $$

where $p_t^{\gamma 1}$, $p_t^{\gamma 2}$ are the transverse momenta of the two photons, $\eta^{\gamma i}$ are their pseudorapidities in the hadronic centre-of-mass frame, and $Y_{\gamma\gamma}$ is the photon-pair rapidity. In the definition of the fiducial volume we do not include the photon-isolation requirement, since this would introduce additional logarithmic corrections of non-global nature in the problem, spoiling the formal N³LL+NNLO accuracy of the differential distributions.<fnref n="4" /> We consider on-shell Higgs boson production followed by a decay into two photons under the narrow-width approximation with a branching ratio of $2.35 \times 10^{-3}$.

In figure 9 we show the comparison of the matched and the fixed-order predictions for the transverse momentum of the photon pair in the fiducial volume, at different perturbative accuracies: N³LL+NLO vs. NLO in the left panel, and N³LL+NNLO vs. NNLO in the right one.

By comparing the two panels of figure 9 we notice a substantial reduction in the theoretical uncertainty in the medium-high-$p_t^{\gamma\gamma}$ region, driven by the increase in perturbative accuracy of the fixed-order computation; at very low $p_t^{\gamma\gamma}$, the prediction is dominated by resummation, which is common to both panels. The pattern observed in the right panel is very similar to what we obtained in the inclusive case in the left panel of figure 7. We

<footnote n="4">⁴However, we point out that photon-isolation criteria in this case are not aggressive, and therefore they could be safely included at fixed order.</footnote>
<figure><img src="image_10.png" /><figcaption>Figure 9. Comparison of the transverse momentum distribution for Higgs boson production at $\sqrt{s}$ = 13 TeV in the fiducial volume defined by eq. (5.1) at N³LL+NLO and NLO (left) and N³LL+NNLO and NNLO (right). The lower panel shows the ratio to the N³LL+NNLO prediction.</figcaption></figure>

stress again that the particularly small uncertainty of the matched prediction is to a certain extent due to the choice of central scales we adopt, namely $\mu_R = \mu_F = m_H/2$, which suffers from large accidental cancellations.

## 6 Results for Drell-Yan production

We now turn to the study of Drell-Yan pair production at the LHC. In this section we present the results for the differential distributions of the transverse momentum of the DY pair, as well as for the angular observable $\phi_{\eta}^{*}$.

We consider 8 TeV proton-proton collisions, and compare the resulting calculation for the differential spectra with ATLAS data from ref. [101]. The fiducial phase-space volume is defined as follows:

$$p_t^{\ell\pm} > 20 \text{ GeV}, \quad |\eta^{\ell\pm}| < 2.4, \quad |Y_{\ell\ell}| < 2.4, \quad 46 \text{ GeV} < M_{\ell\ell} < 150 \text{ GeV}, \quad (6.1)$$

where $p_t^{\ell\pm}$ are the transverse momenta of the two leptons, $\eta^{\ell\pm}$ are their pseudo-rapidities, while $Y_{\ell\ell}$ and $M_{\ell\ell}$ are the rapidity and invariant mass of the di-lepton system, respectively. All rapidities and pseudo-rapidities are evaluated in the hadronic centre-of-mass frame.

For our results, we use parton densities as obtained from the NNPDF3.0 set. The reference value we set for the parameter $p$ appearing in the modified logarithms is $p = 4$, but we have checked that a variation of $p$ by one unit does not induce any significant differences.

We set the central scales as $\mu_R = \mu_F = M_T = \sqrt{M_{\ell\ell}^2 + (p_t^Z)^2}$, while the central resummation scale is chosen to be $x_Q = Q/M_{\ell\ell} = 1/2$. The theoretical uncertainty is estimated through the same set of variations as for Higgs boson production.

The results for $p_t^Z$ and $\phi_{\eta}^*$ are shown in the following two subsections. All plots have the same pattern: the main panels display the comparison of normalised differential distribu-
<figure><img src="image_11.png" /><figcaption>Figure 10. Comparison of the normalised transverse momentum distribution for Drell-Yan pair production at NNLO (green), NNLL+NLO (blue) and N³LL+NNLO (red) at $\sqrt{s}$ = 8 TeV integrated over the full lepton-pair rapidity range (0 < $|Y_{\ell\ell}$| < 2.4), in three different lepton-pair invariant-mass windows. For reference, the ATLAS data is also shown, and the lower panel shows the ratio of each prediction to data.</figcaption></figure>

tions at NNLO (green), NNLL+NLO (blue), and N³LL+NNLO (red), respectively, overlaid on ATLAS data points (black). Correspondingly, the lower insets of each panel show the ratio of the theoretical curves to data, with the same colour code as in the main panels.

## 6.1 Matched predictions for fiducial $p_t^Z$ distributions

In figure 10 we display the normalised $p_t^Z$ distributions in which, in addition to the fiducial cuts reported above, we consider three different lepton-pair invariant-mass windows:

$$ 
\begin{array}{ll}
\text{low invariant mass :} & 46 \text{ GeV} < M_{\ell\ell} < 66 \text{ GeV}, \\
\text{medium invariant mass :} & 66 \text{ GeV} < M_{\ell\ell} < 116 \text{ GeV}, \\
\text{high invariant mass :} & 116 \text{ GeV} < M_{\ell\ell} < 150 \text{ GeV}. \qquad (6.2)
\end{array}
$$
A comparison of the most accurate matched prediction with the fixed-order one shows that the N³LL+NNLO prediction starts differing significantly from the NNLO for $p_t^Z \lesssim 15$ GeV, while for $p_t^Z > 20$ GeV the NNLO is sufficient to provide a reliable description. Comparing matched predictions with different formal accuracy, we note that the N³LL+NNLO curve has a significantly reduced theoretical systematics with respect to that for NNLL+NLO, in the whole $p_t^Z$ range and for all considered invariant-mass windows. The perturbative error is reduced by more than a factor of two at very low $p_t^Z$, where the prediction is dominated by resummation, and the leftover uncertainty in that region is as small as 3-5%, and almost comparable with the excellent experimental precision. The shape of the $p_t^Z$ distributions is also significantly distorted by the inclusion of higher orders: the spectrum is harder than the NNLL+NLO result for $p_t^Z \gtrsim 10$ GeV, and the peak is lower, with the N³LL+NNLO curves in much better agreement with data with respect to NNLL+NLO in the whole kinematic range. Among the three considered windows, the most accurately described at N³LL+NNLO are the ones at intermediate and high invariant mass; the accuracy very slightly degrades for smaller invariant masses, however the theory uncertainty never gets larger than 5-7% over the whole displayed $p_t^Z$ range.

In figure 11 we focus our analysis on the central lepton-pair invariant-mass window defined in eq. (6.2) and show predictions for the normalised $p_t^Z$ distribution in six different lepton-pair rapidity slices:

(a) $0.0 < |Y_{\ell\ell}| < 0.4,$

(b) $0.4 < |Y_{\ell\ell}| < 0.8,$

(c) $0.8 < |Y_{\ell\ell}| < 1.2,$

(d) $1.2 < |Y_{\ell\ell}| < 1.6,$

(e) $1.6 < |Y_{\ell\ell}| < 2.0,$

(f) $2.0 < |Y_{\ell\ell}| < 2.4.$ (6.3)

The comments relevant to figure 10 by far and large apply in this case as well, with our best prediction at N³LL+NNLO affected by an uncertainty that is of order 3-5% in the whole $p_t^Z$ range, regardless of the considered rapidity slice. It is moreover in very good agreement with the experimental data, hence significantly improving on both the NNLL+NLO, in the whole $p_t^Z$ range, and the pure NNLO, in the $p_t^Z \lesssim 20$ GeV region.

## 6.2 Matched predictions for fiducial $\phi_\eta^*$ distributions

Figure 12 shows the $\phi_\eta^*$ distribution for three different lepton-pair invariant-mass windows as defined in eq. (6.2).

The pattern of comparisons among theoretical predictions is qualitatively similar to what discussed for the $p_t^Z$ distribution. Resummation effects at N³LL+NNLO start being important with respect to the pure NNLO in the region $\phi_\eta^* \lesssim 0.2$; the shape of the N³LL+NNLO distribution is significantly distorted with respect to the NNLL+NLO one in a similar fashion as for the $p_t^Z$ case, and the uncertainty band is reduced by a factor of two or more over the whole range and for all invariant-mass windows, down to the level of 3-5% (except at low invariant mass, where the uncertainty is 5-7%).

At variance with the $p_t^Z$ case, however, for $\phi_\eta^*$ we note that the N³LL+NNLO prediction describes data appropriately only in the central- and high-invariant-mass windows. In the low-invariant-mass one, the prediction undershoots data in the medium-hard region, by up to 5-7%. This tension was already observed in the fixed-order NNLO comparison [47].
<figure><img src="image_12.png" /><figcaption>Figure 11. Comparison of the normalised transverse momentum distribution for Drell-Yan pair production at NNLO (green), NNLL+NLO (blue) and N³LL+NNLO (red) at $\sqrt{s} = 8$ TeV in the central lepton-pair invariant-mass window (66 GeV < $M_{\ell\ell}$ < 116 GeV) for six different lepton-pair rapidity slices. For reference, the ATLAS data is also shown, and the lower panel shows the ratio of each prediction to data.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Comparison of the normalised $\phi_{\eta}^{*}$ distribution for Drell-Yan pair production at NNLO (green), NNLL+NLO (blue) and N$^{3}$LL+NNLO (red) at $\sqrt{s} = 8$ TeV integrated over the full lepton-pair rapidity range ($0 < |Y_{\ell\ell}| < 2.4$), in three different lepton-pair invariant-mass windows. For reference, the ATLAS data is also shown, and the lower panel shows the ratio of each prediction to data.</figcaption></figure>

However, given the large statistical uncertainty of the data in this invariant-mass range, the theory still provides a reasonable description of the measurement, and the N$^{3}$LL+NNLO prediction is in much better agreement with data than the NNLL+NLO in the whole range of $\phi_{\eta}^{*}$, especially at low $\phi_{\eta}^{*}$.

In figure 13 we show the results for the $\phi_{\eta}^{*}$ distributions in the central invariant-mass window, see eq. (6.2), split into the six lepton-pair rapidity slices described in eq. (6.3). Moreover, given the availability of experimental measurements, in figures 14 and 15 we also provide predictions sliced in $Y_{\ell\ell}$ for the low- and high- di-lepton invariant-mass windows, respectively. The three rapidity slices we focus on correspond to regions (a+b), (c+d), and (e+f) of eq. (6.3).
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Comparison of the normalised $\phi_{\eta}^*$ distribution for Drell-Yan pair production at NNLO (green), NNLL+NLO (blue) and N$^3$LL+NNLO (red) at $\sqrt{s} = 8$ TeV in the central lepton-pair invariant-mass window (66 GeV $< M_{\ell\ell} < 116$ GeV) for three different lepton-pair rapidity slices. For reference, the ATLAS data is also shown, and the lower panel shows the ratio of each prediction to data.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Comparison of the normalised $\phi_{\eta}^{*}$ distribution for Drell-Yan pair production at NNLO (green), NNLL+NLO (blue) and N$^{3}$LL+NNLO (red) at $\sqrt{s} = 8$ TeV in the low lepton-pair invariant-mass window (46 GeV $< M_{\ell\ell} < 66$ GeV) for three different lepton-pair rapidity slices. For reference, the ATLAS data is also shown, and the lower panel shows the ratio of each prediction to data.</figcaption></figure>

The prediction subdivided in rapidity slices largely shares the same features as that integrated over rapidity, which has been detailed in figure 12. In the central invariant-mass window, data is accurately reproduced by the N$^{3}$LL+NNLO prediction, regardless of the considered rapidity slice, with a theoretical systematics in the 5% range or smaller. The quality of the description slightly degrades at low invariant mass, and to a lesser extent also at high invariant mass, mainly in the hard region, with a pattern similar to that displayed by the rapidity-integrated spectrum. Overall, the uncertainty associated with the N$^{3}$LL+NNLO is of order of 5% or better, with a significant improvement both in the shape and in the systematics with respect to NNLL+NLO.
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Comparison of the normalised $\phi_{\eta}^{*}$ distribution for Drell-Yan pair production at NNLO (green), NNLL+NLO (blue) and N$^{3}$LL+NNLO (red) at $\sqrt{s} = 8$ TeV in the high lepton-pair invariant-mass window ($116 \text{ GeV} < M_{\ell\ell} < 150 \text{ GeV}$) for three different lepton-pair rapidity slices. For reference, the ATLAS data is also shown, and the lower panel shows the ratio of each prediction to data.</figcaption></figure>

# 7 Conclusions

In this work we have presented precise predictions for differential distributions in Higgs boson and Drell-Yan pair production at the LHC at N$^{3}$LL+NNLO.

The resummation is performed in momentum space and is fully exclusive in the Born phase space. For the matching to NNLO we adopted a multiplicative scheme, which allows for the inclusion of the N$^{3}$LO constant terms to the cumulative cross section. These are currently unknown analytically, but can be included numerically once the total N$^{3}$LO cross section has been obtained. The uncertainty associated with the choice of the matching scheme was estimated at NLO accuracy, for which an additive matching with the same ingredients can be also performed. At this order the predictions obtained with the two prescriptions are in very good agreement with each other, and the matching-scheme uncertainty is under control within the perturbative error.
For Higgs boson production in gluon fusion, we have considered the transverse-momentum spectrum both at the inclusive level and in the $H \to \gamma\gamma$ channel within ATLAS fiducial cuts. In both cases, we observe that the resummation reduces the theoretical uncertainties and stabilises the fixed-order result below $p_t^H \sim 40$ GeV. The effects of the N$^3$LL corrections with respect to NNLL+NNLO distributions are moderate in size, with a percent-level correction in most of the range, growing up to 5% at very small $p_t^H$. However, the perturbative uncertainty is reduced significantly below 10 GeV with respect to the NNLL+NNLO case.

For Drell-Yan pair production, we have presented resummed predictions within ATLAS fiducial cuts [101] both for the normalised $p_t^Z$ distributions and for the normalised $\phi_\eta^*$ distributions, and we have compared them to experimental data. In the case of transverse-momentum distributions, the difference between the fixed-order and the N$^3$LL+NNLO result becomes significant for $p_t^Z < 10-15$ GeV, while for $p_t^Z > 20$ GeV the NNLO prediction is sufficient to provide a reliable description of the experimental data. Comparing matched results with different formal accuracy, we note that the N$^3$LL+NNLO prediction has a significantly reduced theoretical uncertainty with respect to that for NNLL+NLO, in the whole $p_t^Z$ range and for all invariant-mass windows considered in our study.

For the $\phi_\eta^*$ distribution, resummation effects start being important with respect to pure NNLO in the region $\phi_\eta^* \lesssim 0.2$. At N$^3$LL+NNLO the shape of the distribution is significantly distorted with respect to that for NNLL+NLO (the spectrum is hardened in the tail, and the height of the peak is lowered), and the uncertainty band is reduced by a factor of two or more over the whole range of $\phi_\eta^*$ and for most invariant-mass windows, down to the level of 3-5%. An exception is at low invariant mass, where the uncertainty remains in the 5-7% range. Unlike the $p_t^Z$ case, for $\phi_\eta^*$ we note that the N$^3$LL+NNLO prediction describes data appropriately only in the central- and high-invariant-mass windows, while at low invariant mass the prediction undershoots the data in the medium-hard region. The difference between the central values of the data and theory here can be of the order of 10%, however no significant tension with the data is observed, due to the sizeable statistical uncertainty in the measurement. The agreement in these invariant-mass bins is much improved by the inclusion of the N$^3$LL+NNLO corrections with respect to the NNLL+NLO distribution.

Our results are an important step in the LHC precision programme, where accurate predictions have become necessary for an appropriate interpretation and exploitation of data. In order to improve on the predictions presented here, several effects must be considered.

For Higgs boson production via gluon fusion, the impact of other heavy quarks, notably the bottom quark, becomes relevant at this level of accuracy and therefore must be taken into account. Recent studies show that the effect of the top-bottom interference at NNLL+NLO [31, 34] could lead to distortions of the transverse-momentum spectrum that are as large as $\sim 5\%$ with respect to the HEFT approximation, and the theory uncertainties associated with this contribution are of $O(20\%)$. These effects are therefore of the same order as the perturbative uncertainties presented here, and must be included for a consistent prediction of the spectrum with 5-10% perturbative accuracy in the region $p_t^H \lesssim m_H$.
In the DY case, the situation is more involved given the smaller perturbative uncertainty. At this level of precision, it is necessary to supplement the predictions obtained in this work at small $p_t^Z$ and $\phi_\eta^*$ with QED corrections and with an estimate of various sources of non-perturbative effects that could be as large as a few % in this region. Similarly, the inclusion of quark masses may have a few-percent effect on the spectrum [147, 148], and more precise studies are necessary in order to assess their impact precisely. Recent analyses [148] suggest that the inclusion of these effects may have a non-negligible impact on observables of current phenomenological interest, such as the determination of the W-boson mass [13]. Given that the size of these effects is of the order of the perturbative uncertainty of the N³LL+NNLO prediction, a careful assessment will be necessary to improve further on the results presented in this work.

## Acknowledgments

XC and TG thank the University of Zurich S3IT for providing the computational resources for this project. XC, TG and AH acknowledge the computing resources provided by the Swiss National Supercomputing Centre (CSCS) under the project ID p501b and UZH10. This research was supported by the Research Executive Agency (REA) of the European Union with the Marie Skłodowska Curie Individual Fellowship contract numbers 702610 (Resummation4PS, PFM) and 659147 (PrecisionTools4LHC, ER) and the ERC Advanced Grant MCatNNLO (340983), by the ERC Consolidator Grant HICCUP (614577), by the ERC Starting Grant PDF4BSM (335260), by the U.K. Science and Technology Facilities Council, and by the Swiss National Science Foundation (SNF) under contracts 200020-175595, 200021-172478 and CRSII2-160814.

## A Formulae for the matching schemes

In this appendix we report the necessary formulae to implement the matching schemes defined in eqs. (4.2) and (4.7) and used in our study. We start by introducing a convenient notation for the perturbative expansion of the various ingredients. We define

$$ \sigma_{\text{tot}}^{\text{N}^3\text{LO}} = \sum_{i=0}^{3} \sigma^{(i)}, \quad \Sigma^{\text{N}^3\text{LO}}(v) = \sigma^{(0)} + \sum_{i=1}^{3} \Sigma^{(i)}(v), \qquad (\text{A.1}) $$

where

$$ \Sigma^{(i)}(v) = \sigma^{(i)} + \bar{\Sigma}^{(i)}(v), \quad \bar{\Sigma}^{(i)}(v) = - \int_{v}^{\infty} dv' \frac{d\Sigma^{(i)}(v')}{dv'}. \qquad (\text{A.2}) $$

Moreover, we denote the perturbative expansion of the resummed cross section $\Sigma^{\text{N}^k\text{LL}}$ as

$$ \Sigma^{\text{EXP}}(v) = \sigma^{(0)} + \sum_{i=1}^{3} \Sigma_{\text{N}^k\text{LL}}^{(i)}(v). \qquad (\text{A.3}) $$
With this notation, the additive scheme of eq. (4.2) becomes (for simplicity we drop the explicit dependence on v in the following)

$$ \Sigma_{\text{add}}^{\text{MAT}} = \Sigma^{\text{N}^k\text{LL}} + \left\{ \sigma^{(1)} + \bar{\Sigma}^{(1)} - \Sigma_{\text{N}^k\text{LL}}^{(1)} \right\} + \left\{ \sigma^{(2)} + \bar{\Sigma}^{(2)} - \Sigma_{\text{N}^k\text{LL}}^{(2)} \right\} + \left\{ \sigma^{(3)} + \bar{\Sigma}^{(3)} - \Sigma_{\text{N}^k\text{LL}}^{(3)} \right\}, \quad (A.4) $$

where the three terms in curly brackets denote the NLO, NNLO and N³LO contributions to the matching, respectively.

For the multiplicative scheme we need to introduce the asymptotic expansion $\Sigma_{\text{asym.}}^{\text{N}^k\text{LL}}$, defined in eq. (4.6) (the definition for $k \neq 3$ is analogous with obvious replacements) in terms of the $\tilde{L} \to 0$ limit of the coefficients $\tilde{\mathcal{L}}_{\text{N}^k\text{LL}}$ of eqs. (3.9), (3.10), (3.11), which read

$$ \tilde{\mathcal{L}}_{\text{NLL}}^{\tilde{L} \to 0} = \sum_{c,c'} \frac{d|M_B|_{cc'}^2}{d\Phi_B} f_c(\mu_F, x_1) f_{c'}(\mu_F, x_2), $$

$$ \tilde{\mathcal{L}}_{\text{NNLL}}^{\tilde{L} \to 0} = \sum_{c,c'} \frac{d|M_B|_{cc'}^2}{d\Phi_B} \sum_{i,j} \int_{x_1}^{1} \frac{dz_1}{z_1} \int_{x_2}^{1} \frac{dz_2}{z_2} f_i \left( \mu_F, \frac{x_1}{z_1} \right) f_j \left( \mu_F, \frac{x_2}{z_2} \right) \\ \times \left\{ \delta_{ci} \delta_{c'j} \delta(1-z_1) \delta(1-z_2) \left( 1 + \frac{\alpha_s(\mu_R)}{2\pi} \tilde{H}^{(1)}(\mu_R, x_Q) \right) \right. \\ \left. + \frac{\alpha_s(\mu_R)}{2\pi} \left( \tilde{C}_{ci}^{(1)}(z_1, \mu_F, x_Q) \delta(1-z_2) \delta_{c'j} + \{z_1 \leftrightarrow z_2; c, i \leftrightarrow c'j\} \right) \right\}, $$

$$ \tilde{\mathcal{L}}_{\text{N}^3\text{LL}}^{\tilde{L} \to 0} = \sum_{c,c'} \frac{d|M_B|_{cc'}^2}{d\Phi_B} \sum_{i,j} \int_{x_1}^{1} \frac{dz_1}{z_1} \int_{x_2}^{1} \frac{dz_2}{z_2} f_i \left( \mu_F, \frac{x_1}{z_1} \right) f_j \left( \mu_F, \frac{x_2}{z_2} \right) \\ \times \left\{ \delta_{ci} \delta_{c'j} \delta(1-z_1) \delta(1-z_2) \left( 1 + \frac{\alpha_s(\mu_R)}{2\pi} \tilde{H}^{(1)}(\mu_R, x_Q) + \frac{\alpha_s^2(\mu_R)}{(2\pi)^2} \tilde{H}^{(2)}(\mu_R, x_Q) \right) \right. \\ + \frac{\alpha_s(\mu_R)}{2\pi} \left( \tilde{C}_{ci}^{(1)}(z_1, \mu_F, x_Q) \delta(1-z_2) \delta_{c'j} + \{z_1 \leftrightarrow z_2; c, i \leftrightarrow c', j\} \right) \\ + \frac{\alpha_s^2(\mu_R)}{(2\pi)^2} \left( \tilde{C}_{ci}^{(2)}(z_1, \mu_F, x_Q) \delta(1-z_2) \delta_{c'j} + \{z_1 \leftrightarrow z_2; c, i \leftrightarrow c', j\} \right) \\ + \frac{\alpha_s^2(\mu_R)}{(2\pi)^2} \left( \tilde{C}_{ci}^{(1)}(z_1, \mu_F, x_Q) \tilde{C}_{c'j}^{(1)}(z_2, \mu_F, x_Q) + G_{ci}^{(1)}(z_1) G_{c'j}^{(1)}(z_2) \right) \\ \left. + \frac{\alpha_s^2(\mu_R)}{(2\pi)^2} \tilde{H}^{(1)}(\mu_R, x_Q) \left( \tilde{C}_{ci}^{(1)}(z_1, \mu_F, x_Q) \delta(1-z_2) \delta_{c'j} + \{z_1 \leftrightarrow z_2; c, i \leftrightarrow c', j\} \right) \right\}. \quad (A.5) $$

In the following formula the perturbative expansion of $\Sigma_{\text{asym.}}^{\text{N}^k\text{LL}}$ is denoted as follows

$$ \Sigma_{\text{asym.}}^{\text{N}^k\text{LL}} = \sigma^{(0)} + \sum_{i=1}^{k-1} \Sigma_{\text{asym.}}^{(i)}. \quad (A.6) $$
With this notation the matching formula (4.7) reads

$$
\begin{equation}
\begin{aligned}
\Sigma_{\text{mult}}^{\text{MAT}}(v) = & \frac{\Sigma_{\text{asym.}}^{\text{N}^k\text{LL}}}{\Sigma_{\text{asym.}}^{\text{N}^k\text{LL}}} \left[ \sigma^{(0)} + \left\{ \sigma^{(1)} + \bar{\Sigma}^{(1)} + \Sigma_{\text{asym.}}^{(1)} - \Sigma_{\text{N}^k\text{LL}}^{(1)} \right\} \right. \\
& \left. \qquad + \left\{ \sigma^{(2)} + \bar{\Sigma}^{(2)} + \Sigma_{\text{asym.}}^{(2)} - \Sigma_{\text{N}^k\text{LL}}^{(2)} + \frac{\Sigma_{\text{asym.}}^{(1)}}{\sigma^{(0)}} (\sigma^{(1)} + \bar{\Sigma}^{(1)}) \right. \right. \\
& \left. \qquad + \frac{(\Sigma_{\text{N}^k\text{LL}}^{(1)})^2}{\sigma^{(0)}} - \frac{\Sigma_{\text{N}^k\text{LL}}^{(1)}}{\sigma^{(0)}} (\sigma^{(1)} + \bar{\Sigma}^{(1)} + \Sigma_{\text{asym.}}^{(1)}) \right\} \\
& + \left\{ \sigma^{(3)} + \bar{\Sigma}^{(3)} - \Sigma_{\text{N}^k\text{LL}}^{(3)} - \frac{(\Sigma_{\text{N}^k\text{LL}}^{(1)})^3}{(\sigma^{(0)})^2} + \frac{(\Sigma_{\text{N}^k\text{LL}}^{(1)})^2}{(\sigma^{(0)})^2} (\sigma^{(1)} + \bar{\Sigma}^{(1)} + \Sigma_{\text{asym.}}^{(1)}) \right. \\
& \left. \qquad + \frac{1}{\sigma_0} \left( (\sigma^{(1)} + \bar{\Sigma}^{(1)})(\Sigma_{\text{asym.}}^{(2)} - \Sigma_{\text{N}^k\text{LL}}^{(2)}) + \Sigma_{\text{asym.}}^{(1)}(\sigma^{(2)} + \bar{\Sigma}^{(2)} - \Sigma_{\text{N}^k\text{LL}}^{(2)}) \right) \right. \\
& \left. \qquad - \frac{1}{(\sigma^{(0)})^2} \Sigma_{\text{N}^k\text{LL}}^{(1)} \left( \Sigma_{\text{asym.}}^{(1)}(\sigma^{(1)} + \bar{\Sigma}^{(1)}) + \sigma^{(0)}(\sigma^{(2)} + \bar{\Sigma}^{(2)} + \Sigma_{\text{asym.}}^{(2)} - 2\Sigma_{\text{N}^k\text{LL}}^{(2)}) \right) \right] ,
\end{aligned}
\tag{A.7}
\end{equation}
$$

where, as above, we grouped the terms entering at NLO, NNLO, and N³LO within curly brackets.

B Formulae for N³LL resummation

In this section we report the expressions for quantities needed for N³LL resummation of
transverse observables, that we have used throughout this article.

First of all we report our convention for the RG equation of the strong coupling
which reads

$$
\frac{d\alpha_s(\mu)}{d\ln\mu^2} = \beta(\alpha_s) \equiv -\alpha_s (\beta_0\alpha_s + \beta_1\alpha_s^2 + \beta_2\alpha_s^3 + \beta_3\alpha_s^4 + \dots), \quad (B.1)
$$

where the coefficients of the β-function are

$$
\beta_0 = \frac{11C_A - 2n_f}{12\pi}, \quad \beta_1 = \frac{17C_A^2 - 5C_A n_f - 3C_F n_f}{24\pi^2}, \qquad (B.2)
$$

$$
\beta_2 = \frac{2857 C_A^3 + (54 C_F^2 - 615 C_F C_A - 1415 C_A^2) n_f + (66 C_F + 79 C_A) n_f^2}{3456 \pi^3}, \quad (B.3)
$$

$$
\begin{equation}
\beta_3 = \frac{1}{(4\pi)^4} \left\{
\begin{aligned}[t]
& C_A C_F n_f^2 \frac{1}{4} \left( \frac{17152}{243} + \frac{448}{9} \zeta_3 \right) + C_A C_F^2 n_f \frac{1}{2} \left( -\frac{4204}{27} + \frac{352}{9} \zeta_3 \right) \\
& + \frac{53}{243} C_A n_f^3 + C_A^2 C_F n_f \frac{1}{2} \left( \frac{7073}{243} - \frac{656}{9} \zeta_3 \right) + C_A^2 n_f^2 \frac{1}{4} \left( \frac{7930}{81} + \frac{224}{9} \zeta_3 \right) \\
& + \frac{154}{243} C_F n_f^3 + C_A^3 n_f \frac{1}{2} \left( -\frac{39143}{81} + \frac{136}{3} \zeta_3 \right) + C_A^4 \left( \frac{150653}{486} - \frac{44}{9} \zeta_3 \right) \\
& + C_F^2 n_f^2 \frac{1}{4} \left( \frac{1352}{27} - \frac{704}{9} \zeta_3 \right) + 23 C_F^3 n_f + n_f \frac{d_F^{abcd} d_A^{abcd}}{N_A} \left( \frac{512}{9} - \frac{1664}{3} \zeta_3 \right) \\
& + n_f^2 \frac{d_F^{abcd} d_F^{abcd}}{N_A} \left( -\frac{704}{9} + \frac{512}{3} \zeta_3 \right) + \frac{d_A^{abcd} d_A^{abcd}}{N_A} \left( -\frac{80}{9} + \frac{704}{3} \zeta_3 \right)
\end{aligned}
\right\}, \tag{B.4}
\end{equation}
$$
with

$$ \frac{d_F^{abcd} d_F^{abcd}}{N_A} = \frac{N_c^4 - 6N_c^2 + 18}{96N_c^2}, \quad \frac{d_F^{abcd} d_A^{abcd}}{N_A} = \frac{N_c(N_c^2 + 6)}{48}, \quad \frac{d_A^{abcd} d_A^{abcd}}{N_A} = \frac{N_c^2(N_c^2 + 36)}{24}, $$

and $C_A = N_c$, $C_F = \frac{N_c^2-1}{2N_c}$, and $N_c = 3$.

We also provide expressions for the functions $g_i(\lambda)$ entering in the N$^3$LL Sudakov radiator eq. (3.6) and its derivative. We define

$$ \lambda = \alpha_s(\mu_R)\beta_0 \tilde{L}. \qquad (B.5) $$

We have:

$$ g_1(\lambda) = \frac{A^{(1)} 2\lambda + \ln(1 - 2\lambda)}{\pi \beta_0 2\lambda}, \qquad (B.6) $$

$$ g_2(\lambda) = \frac{1}{2\pi\beta_0} \ln(1-2\lambda) \left( A^{(1)} \ln \frac{1}{x_Q^2} + B^{(1)} \right) - \frac{A^{(2)}}{4\pi^2\beta_0^2} \frac{2\lambda + (1-2\lambda)\ln(1-2\lambda)}{1-2\lambda} \\ + A^{(1)} \left( -\frac{\beta_1 \ln(1-2\lambda)((2\lambda-1)\ln(1-2\lambda)-2) - 4\lambda}{1-2\lambda} \right. \\ \left. - \frac{1}{2\pi\beta_0} \frac{(2\lambda(1-\ln(1-2\lambda)) + \ln(1-2\lambda))}{1-2\lambda} \ln \frac{\mu_R^2}{x_Q^2 M^2} \right), \qquad (B.7) $$

$$ g_3(\lambda) = \left( A^{(1)} \ln \frac{1}{x_Q^2} + B^{(1)} \right) \left( -\frac{\lambda}{1-2\lambda} \ln \frac{\mu_R^2}{x_Q^2 M^2} + \frac{\beta_1}{2\beta_0^2} \frac{2\lambda + \ln(1-2\lambda)}{1-2\lambda} \right) \\ - \frac{1}{2\pi\beta_0} \frac{\lambda}{1-2\lambda} \left( A^{(2)} \ln \frac{1}{x_Q^2} + B^{(2)} \right) - \frac{A^{(3)}}{4\pi^2\beta_0^2} \frac{\lambda^2}{(1-2\lambda)^2} \\ + A^{(2)} \left( \frac{\beta_1}{4\pi\beta_0^3} \frac{2\lambda(3\lambda-1) + (4\lambda-1)\ln(1-2\lambda)}{(1-2\lambda)^2} - \frac{1}{\pi\beta_0} \frac{\lambda^2}{(1-2\lambda)^2} \ln \frac{\mu_R^2}{x_Q^2 M^2} \right) \\ + A^{(1)} \left( \frac{\lambda(\beta_0\beta_2(1-3\lambda) + \beta_1^2\lambda)}{\beta_0^4(1-2\lambda)^2} + \frac{(1-2\lambda)\ln(1-2\lambda)(\beta_0\beta_2(1-2\lambda) + 2\beta_1^2\lambda)}{2\beta_0^4(1-2\lambda)^2} \right. \\ \left. + \frac{\beta_1^2}{4\beta_0^4} \frac{(1-4\lambda)\ln^2(1-2\lambda)}{(1-2\lambda)^2} - \frac{\lambda^2}{(1-2\lambda)^2} \ln^2 \frac{\mu_R^2}{x_Q^2 M^2} \right. \\ \left. - \frac{\beta_1}{2\beta_0^2} \frac{(2\lambda(1-2\lambda) + (1-4\lambda)\ln(1-2\lambda))}{(1-2\lambda)^2} \ln \frac{\mu_R^2}{x_Q^2 M^2} \right), \qquad (B.8) $$

$$ g_4(\lambda) = \frac{A^{(4)}(3-2\lambda)\lambda^2}{24\pi^2\beta_0^2(2\lambda-1)^3} \\ + \frac{A^{(3)}}{48\pi\beta_0^3(2\lambda-1)^3} \left\{ 3\beta_1(1-6\lambda)\ln(1-2\lambda) + 2\lambda \left( \beta_1(5\lambda(2\lambda-3)+3) \right. \right. \\ \left. \left. + 6\beta_0^2(3-2\lambda)\lambda \ln \frac{\mu_R^2}{x_Q^2 M^2} \right) + 12\beta_0^2(\lambda-1)\lambda(2\lambda-1) \ln \frac{1}{x_Q^2} \right\} $$
$$
\begin{equation}
\begin{aligned}
& + \frac{A^{(2)}}{24\beta_0^4(2\lambda-1)^3} \left\{
\begin{aligned}[t]
& 32\beta_0\beta_2\lambda^3 - 2\beta_1^2\lambda(\lambda(22\lambda-9) + 3) \\
& + 12\beta_0^4(3-2\lambda)\lambda^2\ln^2\frac{\mu_R^2}{x_Q^2 M^2} + 6\beta_0^2\ln\frac{\mu_R^2}{x_Q^2 M^2} \times \\
& \quad \left( \beta_1(1-6\lambda)\ln(1-2\lambda) + 2(\lambda-1)\lambda(2\lambda-1) \right) \left( \beta_1 + 2\beta_0^2\ln\frac{1}{x_Q^2} \right) \\
& + 3\beta_1 \left( \beta_1\ln(1-2\lambda)(2\lambda + (6\lambda-1)\ln(1-2\lambda) - 1) \right. \\
& \qquad \left. - 2\beta_0^2(2\lambda-1)(2(\lambda-1)\lambda - \ln(1-2\lambda))\ln\frac{1}{x_Q^2} \right)
\end{aligned}
\right\} \\
& + \frac{\pi A^{(1)}}{12\beta_0^5(2\lambda-1)^3} \left\{
\begin{aligned}[t]
& \beta_1^3(1-6\lambda)\ln^3(1-2\lambda) + 3\ln(1-2\lambda) \left( \beta_0^2\beta_3(2\lambda-1)^3 \right. \\
& \qquad \left. + \beta_0\beta_1\beta_2(1-2\lambda(8\lambda^2-4\lambda+3)) + 4\beta_1^3\lambda^2(2\lambda+1) \right) \\
& + \beta_0^2\beta_1\ln\frac{\mu_R^2}{x_Q^2 M^2} \left( \beta_0^2(1-6\lambda)\ln\frac{\mu_R^2}{x_Q^2 M^2} - 4\beta_1\lambda \right) \\
& + 3\beta_1^2\ln^2(1-2\lambda) \left( 2\beta_1\lambda + \beta_0^2(6\lambda-1)\ln\frac{\mu_R^2}{x_Q^2 M^2} \right) \\
& + 3\beta_0^2(2\lambda-1)\ln\frac{1}{x_Q^2} \left( -\beta_1^2\ln^2(1-2\lambda) + 2\beta_0^2\beta_1\ln(1-2\lambda)\ln\frac{\mu_R^2}{x_Q^2 M^2} \right. \\
& \qquad \left. + 4\lambda \left( \lambda(\beta_1^2 - \beta_0\beta_2) + \beta_0^4(\lambda-1)\ln^2\frac{\mu_R^2}{x_Q^2 M^2} \right) \right) \\
& + 2\lambda \left( \beta_0^2\beta_3((15-14\lambda)\lambda-3) + \beta_0\beta_1\beta_2(5\lambda(2\lambda-3)+3) \right. \\
& \qquad \left. + 4\beta_1^3\lambda^2 + 2\beta_0^6(3-2\lambda)\lambda\ln^3\frac{\mu_R^2}{x_Q^2 M^2} + 3\beta_0^4\beta_1\ln^2\frac{\mu_R^2}{x_Q^2 M^2} \right. \\
& \qquad \left. + 6\beta_0^2\lambda(2\lambda+1)(\beta_0\beta_2 - \beta_1^2)\ln\frac{\mu_R^2}{x_Q^2 M^2} - 8\beta_0^6(4\lambda^2-6\lambda+3)\zeta_3 \right)
\end{aligned}
\right\} \\
& + \frac{B^{(3)}(\lambda-1)\lambda}{4\pi\beta_0(1-2\lambda)^2} + \frac{B^{(2)}\left(\beta_1\ln(1-2\lambda) - 2(\lambda-1)\lambda\left(\beta_1 - 2\beta_0^2\ln\frac{\mu_R^2}{x_Q^2 M^2}\right)\right)}{4\beta_0^2(1-2\lambda)^2} \\
& + \frac{\pi B^{(1)}}{4\beta_0^3(1-2\lambda)^2} \left\{
\begin{aligned}[t]
& 4\lambda \left( \lambda(\beta_1^2 - \beta_0\beta_2) + \beta_0^4(\lambda-1)\ln^2\frac{\mu_R^2}{x_Q^2 M^2} \right) \\
& - \beta_1^2\ln^2(1-2\lambda) + 2\beta_0^2\beta_1\ln(1-2\lambda)\ln\frac{\mu_R^2}{x_Q^2 M^2}
\end{aligned}
\right\}.
\end{aligned}
\tag{B.9}
\end{equation}
$$
For Higgs boson production in gluon fusion, the coefficients $A^{(i)}$ and $B^{(i)}$ which enter the formulae above are (in units of $\alpha_s/(2\pi)$)

$$A_{\text{ggH}}^{(1)} = 2C_A,$$

$$A_{\text{ggH}}^{(2)} = \left(\frac{67}{9} - \frac{\pi^2}{3}\right) C_A^2 - \frac{10}{9} C_A n_f,$$

$$A_{\text{ggH}}^{(3)} = \left(-22\zeta_3 - \frac{67\pi^2}{27} + \frac{11\pi^4}{90} + \frac{15503}{324}\right) C_A^3 + \left(\frac{10\pi^2}{27} - \frac{2051}{162}\right) C_A^2 n_f \\ + \left(4\zeta_3 - \frac{55}{12}\right) C_A C_F n_f + \frac{50}{81} C_A n_f^2,$$

$$A_{\text{ggH}}^{(4)} = \left(\frac{121}{3}\zeta_3\zeta_2 - \frac{8789\zeta_2}{162} - \frac{19093\zeta_3}{54} - \frac{847\zeta_4}{24} + 132\zeta_5 + \frac{3761815}{11664}\right) C_A^4 + \left(-\frac{4\zeta_3}{9} - \frac{232}{729}\right) C_A n_f^3 \\ + \left(-\frac{22}{3}\zeta_3\zeta_2 + \frac{2731\zeta_2}{162} + \frac{4955\zeta_3}{54} + \frac{11\zeta_4}{6} - 24\zeta_5 - \frac{31186}{243}\right) C_A^3 n_f \\ + \left(-\frac{38\zeta_3}{9} - 2\zeta_4 + \frac{215}{24}\right) C_A C_F n_f^2 + \left(\frac{272\zeta_3}{9} + 11\zeta_4 - \frac{7351}{144}\right) C_A^2 C_F n_f \\ + \left(-\frac{103\zeta_2}{81} - \frac{47\zeta_3}{27} + \frac{5\zeta_4}{6} + \frac{13819}{972}\right) C_A^2 n_f^2 + \Gamma_{\text{cusp,ggH}}^{(4)} + C_A \Delta A^{(4)},$$

$$B_{\text{ggH}}^{(1)} = -\frac{11}{3}C_A + \frac{2}{3}n_f,$$

$$B_{\text{ggH}}^{(2)} = \left(\frac{11\zeta_2}{6} - 6\zeta_3 - \frac{16}{3}\right) C_A^2 + \left(\frac{4}{3} - \frac{\zeta_2}{3}\right) C_A n_f + n_f C_F,$$

$$B_{\text{ggH}}^{(3)} = \left(\frac{22\zeta_3\zeta_2}{3} - \frac{799\zeta_2}{81} - \frac{5\pi^2\zeta_3}{9} - \frac{2533\zeta_3}{54} - \frac{77\zeta_4}{12} + 20\zeta_5 - \frac{319\pi^4}{1080} + \frac{6109\pi^2}{1944} + \frac{34219}{1944}\right) C_A^3 \\ + \left(\frac{103\zeta_2}{81} + \frac{202\zeta_3}{27} - \frac{5\zeta_4}{6} + \frac{41\pi^4}{540} - \frac{599\pi^2}{972} - \frac{10637}{1944}\right) C_A^2 n_f \\ + \left(-\frac{2\zeta_3}{27} + \frac{5\pi^2}{162} + \frac{529}{1944}\right) C_A n_f^2 + \left(2\zeta_4 - \frac{\pi^4}{45} - \frac{\pi^2}{12} + \frac{241}{72}\right) C_A C_F n_f \\ - \frac{1}{4} C_F^2 n_f - \frac{11}{36} C_A n_f^2 + C_A \Delta B^{(3)}. \qquad (B.10)$$

For Drell-Yan production, the coefficients read

$$A_{\text{DY}}^{(1)} = 2C_F,$$

$$A_{\text{DY}}^{(2)} = \left(\frac{67}{9} - \frac{\pi^2}{3}\right) C_A C_F - \frac{10}{9} C_F n_f,$$

$$A_{\text{DY}}^{(3)} = \left(\frac{15503}{324} - \frac{67\pi^2}{27} + \frac{11\pi^4}{90} - 22\zeta_3\right) C_A^2 C_F + \left(-\frac{2051}{162} + \frac{10\pi^2}{27}\right) C_A C_F n_f \\ + \left(-\frac{55}{12} + 4\zeta_3\right) C_F^2 n_f + \frac{50}{81} C_F n_f^2,$$

$$A_{\text{DY}}^{(4)} = \left(\frac{3761815}{11664} - \frac{8789\zeta_2}{162} - \frac{19093\zeta_3}{54} + \frac{121\zeta_2\zeta_3}{3} - \frac{847\zeta_4}{24} + 132\zeta_5\right) C_A^3 C_F \\ + \left(-\frac{232}{729} - \frac{4\zeta_3}{9}\right) C_F n_f^3 + \left(\frac{215}{24} - \frac{38\zeta_3}{9} - 2\zeta_4\right) C_F^2 n_f^2$$
$$ 
\begin{aligned}
& + \left( -\frac{31186}{243} + \frac{2731\zeta_2}{162} + \frac{4955\zeta_3}{54} - \frac{22\zeta_2\zeta_3}{3} + \frac{11\zeta_4}{6} - 24\zeta_5 \right) C_A^2 C_F n_f \\
& + \left( -\frac{7351}{144} + \frac{272\zeta_3}{9} + 11\zeta_4 \right) C_A C_F^2 n_f + \left( \frac{13819}{972} - \frac{103\zeta_2}{81} - \frac{47\zeta_3}{27} + \frac{5\zeta_4}{4} \right) C_A C_F n_f^2 \\
& + \Gamma_{\text{cusp,DY}}^{(4)} + C_F \Delta A^{(4)},
\end{aligned}
$$

$$B_{\text{DY}}^{(1)} = -3C_F,$$

$$B_{\text{DY}}^{(2)} = \left( -\frac{17}{12} - \frac{11\pi^2}{12} + 6\zeta_3 \right) C_A C_F + \left( -\frac{3}{4} + \pi^2 - 12\zeta_3 \right) C_F^2 + \left( \frac{1}{6} + \frac{\pi^2}{6} \right) C_F n_f,$$

$$ 
\begin{aligned}
B_{\text{DY}}^{(3)} = & \left( \frac{22\zeta_3\zeta_2}{3} - \frac{799\zeta_2}{81} - \frac{11\pi^2\zeta_3}{9} + \frac{2207\zeta_3}{54} - \frac{77\zeta_4}{12} - 10\zeta_5 - \frac{83\pi^4}{360} - \frac{7163\pi^2}{1944} + \frac{151571}{3888} \right) C_A^2 C_F \\
& + \left( \frac{4\pi^2 - 51}{3}\zeta_3 + 60\zeta_5 - \frac{2\pi^4}{5} - \frac{3\pi^2}{4} - \frac{29}{8} \right) C_F^3 + \left( \frac{34\zeta_3}{3} + 2\zeta_4 - \frac{7\pi^4}{54} - \frac{13\pi^2}{36} + \frac{23}{4} \right) C_F^2 n_f \\
& + \left( -\frac{2}{3}\pi^2\zeta_3 - \frac{211\zeta_3}{3} - 30\zeta_5 + \frac{247\pi^4}{540} + \frac{205\pi^2}{36} - \frac{151}{16} \right) C_A C_F^2 \\
& + \left( \frac{103\zeta_2}{81} - \frac{128\zeta_3}{27} - \frac{5\zeta_4}{6} + \frac{11\pi^4}{180} + \frac{1297\pi^2}{972} - \frac{3331}{243} \right) C_A C_F n_f \\
& + \left( \frac{10\zeta_3}{27} - \frac{5\pi^2}{54} + \frac{1115}{972} \right) C_F n_f^2 + C_F \Delta B^{(3)}.
\end{aligned}
\quad (B.11) $$

The expressions for the coefficients $A^{(i)}$ and $B^{(i)}$ are extracted from refs. [62, 92, 93, 149] for Higgs boson production and refs. [67, 92, 93, 150] for DY production. The $N^3\text{LL}$ anomalous dimension $A^{(4)}$ receives a contribution from the four-loop cusp anomalous dimension $\Gamma_{\text{cusp}}^{(4)}$, that has recently been computed numerically in ref. [95], and is given by

$$ 
\begin{aligned}
\Gamma_{\text{cusp,ggH}}^{(4)} & \approx 2555 - 732.125 n_f + 27.5031 n_f^2 + 0.460173 n_f^3, \\
\Gamma_{\text{cusp,DY}}^{(4)} & \approx 1293.88 - 323.244 n_f + 12.2236 n_f^2 + 0.204522 n_f^3.
\end{aligned}
\quad (B.12) $$

The extra terms

$$ \Delta A^{(4)} = -64\pi^3\beta_0^3\zeta_3, \quad \Delta B^{(3)} = -32\pi^2\beta_0^2\zeta_3, \quad \Delta H^{(2)} = \frac{16}{3}\pi\beta_0\zeta_3, \quad (B.13) $$

are a feature of performing the resummation in momentum space, and do not appear in the anomalous dimensions in b space (see ref. [85] for details). The term $\Delta H^{(2)}$ will appear in the $\tilde{H}$ functions defined below.

We also present the expansion of hard-virtual coefficient function $H$ in powers of the strong coupling

$$ H(M) = 1 + \sum_{n=1}^{2} \left( \frac{\alpha_s(M)}{2\pi} \right)^n H^{(n)}(M), \quad (B.14) $$

with

$$ 
\begin{aligned}
H_{\text{ggH}}^{(1)}(M) &= C_A \left( 5 + \frac{7}{6}\pi^2 \right) - 3C_F, \\
H_{\text{ggH}}^{(2)}(M) &= \frac{5359}{54} + \frac{137}{6} \ln \frac{m_H^2}{m_t^2} + \frac{1679}{24}\pi^2 + \frac{37}{8}\pi^4 - \frac{499}{6}\zeta_3 + C_A \Delta H^{(2)}, \quad n_f = 5, \quad (B.15)
\end{aligned}
$$
and

$$ 
\begin{aligned}
H_{\text{DY}}^{(1)}(M) &= C_F \left(5 + \frac{7}{6}\pi^2\right), \\
H_{\text{DY}}^{(2)}(M) &= -\frac{57433}{972} + \frac{281}{162}\pi^2 + \frac{22}{27}\pi^4 + \frac{1178}{27}\zeta_3 + C_F \Delta H^{(2)}, \quad n_f = 5.
\end{aligned} 
\quad (B.16) 
$$

The factors $\tilde{H}$ that appear in the luminosity prefactors (eqs. (3.9), (3.10), (3.11)) are defined as

$$ 
\begin{aligned}
\tilde{H}^{(1)}(\mu_R, x_Q) &= H^{(1)}(\mu_R) + \left(-\frac{1}{2}A^{(1)}\ln x_Q^2 + B^{(1)}\right)\ln x_Q^2, \\
\tilde{H}^{(2)}(\mu_R, x_Q) &= H^{(2)}(\mu_R) + \frac{(A^{(1)})^2}{8}\ln^4 x_Q^2 - \left(\frac{A^{(1)}B^{(1)}}{2} + \frac{A^{(1)}}{3}\pi\beta_0\right)\ln^3 x_Q^2 \\
&\quad + \left(\frac{-A^{(2)} + (B^{(1)})^2}{2} + \pi\beta_0\left(B^{(1)} + A^{(1)}\ln\frac{x_Q^2 M^2}{\mu_R^2}\right)\right)\ln^2 x_Q^2 \\
&\quad - \left(-B^{(2)} + B^{(1)}2\pi\beta_0\ln\frac{x_Q^2 M^2}{\mu_R^2}\right)\ln x_Q^2 + H^{(1)}(\mu_R)\ln x_Q^2\left(-\frac{1}{2}A^{(1)}\ln x_Q^2 + B^{(1)}\right).
\end{aligned} 
\quad (B.17) 
$$

Finally we report the expansion of the collinear coefficient functions $C_{ab}$

$$ C_{ab}(z) = \delta(1-z)\delta_{ab} + \sum_{n=1}^{2} \left(\frac{\alpha_s(\mu)}{2\pi}\right)^n C_{ab}^{(n)}(z), \quad (B.18) $$

where $\mu$ is the same scale that enters parton densities. The first-order expansion has been known for a long time and reads

$$ C_{ab}^{(1)}(z) = -\hat{P}_{ab}^{(0),\epsilon}(z) - \delta_{ab}\delta(1-z)\frac{\pi^2}{12}, \quad (B.19) $$

where $\hat{P}_{ab}^{(0),\epsilon}(z)$ is the $O(\epsilon)$ part of the leading-order regularised splitting functions $\hat{P}_{ab}^{(0)}(z)$

$$ 
\begin{aligned}
\hat{P}_{qq}^{(0)}(z) &= C_F \left[ \frac{1+z^2}{(1-z)_+} + \frac{3}{2}\delta(1-z) \right], & \hat{P}_{qq}^{(0),\epsilon}(z) &= -C_F(1-z), \\
\hat{P}_{qg}^{(0)}(z) &= \frac{1}{2}[z^2+(1-z)^2], & \hat{P}_{qg}^{(0),\epsilon}(z) &= -z(1-z), \\
\hat{P}_{gq}^{(0)}(z) &= C_F \frac{1+(1-z)^2}{z}, & \hat{P}_{gq}^{(0),\epsilon}(z) &= -C_F z, \\
\hat{P}_{gg}^{(0)}(z) &= 2C_A \left[ \frac{z}{(1-z)_+} + \frac{1-z}{z} + z(1-z) \right] + 2\pi\beta_0\delta(1-z), & \hat{P}_{gg}^{(0),\epsilon}(z) &= 0.
\end{aligned} 
\quad (B.20) 
$$

The second-order collinear coefficient functions $C_{ab}^{(2)}(z)$, as well as the $G$ coefficients (see eqs. (3.9), (3.10), (3.11)) for gluon-fusion processes are obtained in refs. [88, 90, 91], while for quark-induced processes they are derived in ref. [89]. In the present work we extract their expressions using the results of refs. [88, 89]. For gluon-fusion processes, the $C_{gq}^{(2)}$
and $C_{gg}^{(2)}$ coefficients normalised as in eq. (B.19) are extracted from eqs. (30) and (32) of ref. [88], respectively, where we use the hard coefficients of eqs. (B.15) without the new term $\Delta H^{(2)}$ in the $H_g^{(2)}(M)$ coefficient.<fnref n="5" /> The coefficient $G^{(1)}$ is taken from eq. (13) of ref. [88]. Similarly, for quark-initiated processes, we extract $C_{qg}^{(2)}$ and $C_{qq}^{(2)}$ from eqs. (32) and (34) of ref. [89], respectively, where we use the hard coefficients from eqs. (B.16) without the new term $\Delta H^{(2)}$ in the $H_q^{(2)}(M)$ coefficient. The remaining quark coefficient function $\tilde{C}_{q\bar{q}}^{(2)}$, $C_{q\bar{q}'}^{(2)}$ and $C_{qq'}^{(2)}$ are extracted from eq. (35) of the same article.

The coefficients $\tilde{C}$ in eqs. (3.9), (3.10), (3.11) are defined as

$$ 
\begin{aligned}
\tilde{C}_{ab}^{(1)}(z, \mu_F, x_Q) &= C_{ab}^{(1)}(z) + \hat{P}_{ab}^{(0)}(z) \ln \frac{x_Q^2 M^2}{\mu_F^2}, \\
\tilde{C}_{ab}^{(2)}(z, \mu_F, x_Q) &= C_{ab}^{(2)}(z) + \pi \beta_0 \hat{P}_{ab}^{(0)}(z) \left( \ln^2 \frac{x_Q^2 M^2}{\mu_F^2} - 2 \ln \frac{x_Q^2 M^2}{\mu_F^2} \ln \frac{x_Q^2 M^2}{\mu_R^2} \right) \\
&\quad + \hat{P}_{ab}^{(1)}(z) \ln \frac{x_Q^2 M^2}{\mu_F^2} + \frac{1}{2} (\hat{P}^{(0)} \otimes \hat{P}^{(0)})_{ab}(z) \ln^2 \frac{x_Q^2 M^2}{\mu_F^2} \\
&\quad + (C^{(1)} \otimes \hat{P}^{(0)})_{ab}(z) \ln \frac{x_Q^2 M^2}{\mu_F^2} - 2 \pi \beta_0 C_{ab}^{(1)}(z) \ln \frac{x_Q^2 M^2}{\mu_R^2}.
\end{aligned}
\quad (B.21) $$

**Open Access.** This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.

## References

[1] CMS collaboration, Measurement of the transverse momentum spectra of weak vector bosons produced in proton-proton collisions at $\sqrt{s} = 8$ TeV, JHEP 02 (2017) 096 [arXiv:1606.05864] [INSPIRE].

[2] LHCb collaboration, Measurement of the forward Z boson production cross-section in pp collisions at $\sqrt{s} = 13$ TeV, JHEP 09 (2016) 136 [arXiv:1607.06495] [INSPIRE].

[3] ATLAS collaboration, Measurement of differential cross sections and $W^+/W^-$ cross-section ratios for W boson production in association with jets at $\sqrt{s} = 8$ TeV with the ATLAS detector, JHEP 05 (2018) 077 [arXiv:1711.03296] [INSPIRE].

[4] ATLAS collaboration, Measurement of inclusive and differential cross sections in the $H \to ZZ^* \to 4\ell$ decay channel in pp collisions at $\sqrt{s} = 13$ TeV with the ATLAS detector, JHEP 10 (2017) 132 [arXiv:1708.02810] [INSPIRE].

[5] CMS collaboration, Measurement of differential cross sections in the kinematic angular variable $\phi^*$ for inclusive Z boson production in pp collisions at $\sqrt{s} = 8$ TeV, JHEP 03 (2018) 172 [arXiv:1710.07955] [INSPIRE].

[6] CMS collaboration, Measurements of properties of the Higgs boson decaying into the four-lepton final state in pp collisions at $\sqrt{s} = 13$ TeV, JHEP 11 (2017) 047 [arXiv:1706.09936] [INSPIRE].

<footnote n="5">⁵These must be replaced by $H^{(1)} \to H^{(1)}/2$ and $H^{(2)} \to H^{(2)}/4$ to match the convention of refs. [88, 89].</footnote>
[7] ATLAS collaboration, *Measurements of Higgs boson properties in the diphoton decay channel with 36 fb⁻¹ of pp collision data at √s = 13 TeV* with the ATLAS detector, Phys. Rev. D 98 (2018) 052005 [arXiv:1802.04146] [INSPIRE].

[8] CMS collaboration, *Measurements of Higgs boson properties in the diphoton decay channel in proton-proton collisions at √s = 13 TeV*, JHEP 11 (2018) 185 [arXiv:1804.02716] [INSPIRE].

[9] NNPDF collaboration, R.D. Ball et al., Parton distributions from high-precision collider data, Eur. Phys. J. C 77 (2017) 663 [arXiv:1706.00428] [INSPIRE].

[10] I. Scimemi and A. Vladimirov, Power corrections and renormalons in Transverse Momentum Distributions, JHEP 03 (2017) 002 [arXiv:1609.06047] [INSPIRE].

[11] I. Scimemi and A. Vladimirov, Analysis of vector boson production within TMD factorization, Eur. Phys. J. C 78 (2018) 89 [arXiv:1706.01473] [INSPIRE].

[12] A. Bacchetta, F. Delcarro, C. Pisano, M. Radici and A. Signori, Extraction of partonic transverse momentum distributions from semi-inclusive deep-inelastic scattering, Drell-Yan and Z-boson production, JHEP 06 (2017) 081 [arXiv:1703.10157] [INSPIRE].

[13] ATLAS collaboration, Measurement of the W-boson mass in pp collisions at √s = 7 TeV with the ATLAS detector, Eur. Phys. J. C 78 (2018) 110 [Erratum ibid. C 78 (2018) 898] [arXiv:1701.07240] [INSPIRE].

[14] R. Boughezal, A. Guffanti, F. Petriello and M. Ubiali, The impact of the LHC Z-boson transverse momentum data on PDF determinations, JHEP 07 (2017) 130 [arXiv:1705.00343] [INSPIRE].

[15] F. Bishara, U. Haisch, P.F. Monni and E. Re, Constraining Light-Quark Yukawa Couplings from Higgs Distributions, Phys. Rev. Lett. 118 (2017) 121801 [arXiv:1606.09253] [INSPIRE].

[16] Y. Soreq, H.X. Zhu and J. Zupan, Light quark Yukawa couplings from Higgs kinematics, JHEP 12 (2016) 045 [arXiv:1606.09621] [INSPIRE].

[17] A. Banfi, S. Redford, M. Vesterinen, P. Waller and T.R. Wyatt, Optimisation of variables for studying dilepton transverse momentum distributions at hadron colliders, Eur. Phys. J. C 71 (2011) 1600 [arXiv:1009.1580] [INSPIRE].

[18] D. de Florian, M. Grazzini and Z. Kunszt, Higgs production with large transverse momentum in hadronic collisions at next-to-leading order, Phys. Rev. Lett. 82 (1999) 5209 [hep-ph/9902483] [INSPIRE].

[19] R.V. Harlander and W.B. Kilgore, Next-to-next-to-leading order Higgs production at hadron colliders, Phys. Rev. Lett. 88 (2002) 201801 [hep-ph/0201206] [INSPIRE].

[20] C. Anastasiou and K. Melnikov, Higgs boson production at hadron colliders in NNLO QCD, Nucl. Phys. B 646 (2002) 220 [hep-ph/0207004] [INSPIRE].

[21] V. Ravindran, J. Smith and W.L. van Neerven, NNLO corrections to the total cross-section for Higgs boson production in hadron hadron collisions, Nucl. Phys. B 665 (2003) 325 [hep-ph/0302135] [INSPIRE].

[22] V. Ravindran, J. Smith and W.L. van Neerven, Next-to-leading order QCD corrections to differential distributions of Higgs boson production in hadron hadron collisions, Nucl. Phys. B 634 (2002) 247 [hep-ph/0201114] [INSPIRE].
[23] C. Anastasiou et al., High precision determination of the gluon fusion Higgs boson cross-section at the LHC, JHEP 05 (2016) 058 [arXiv:1602.00695] [INSPIRE].

[24] B. Mistlberger, Higgs boson production at hadron colliders at N³LO in QCD, JHEP 05 (2018) 028 [arXiv:1802.00833] [INSPIRE].

[25] M. Bonvini and S. Marzani, Double resummation for Higgs production, Phys. Rev. Lett. 120 (2018) 202003 [arXiv:1802.07758] [INSPIRE].

[26] R. Boughezal, F. Caola, K. Melnikov, F. Petriello and M. Schulze, Higgs boson production in association with a jet at next-to-next-to-leading order, Phys. Rev. Lett. 115 (2015) 082003 [arXiv:1504.07922] [INSPIRE].

[27] R. Boughezal, C. Focke, W. Giele, X. Liu and F. Petriello, Higgs boson production in association with a jet at NNLO using jettiness subtraction, Phys. Lett. B 748 (2015) 5 [arXiv:1505.03893] [INSPIRE].

[28] F. Caola, K. Melnikov and M. Schulze, Fiducial cross sections for Higgs boson production in association with a jet at next-to-next-to-leading order in QCD, Phys. Rev. D 92 (2015) 074032 [arXiv:1508.02684] [INSPIRE].

[29] X. Chen, J. Cruz-Martinez, T. Gehrmann, E.W.N. Glover and M. Jaquier, NNLO QCD corrections to Higgs boson production at large transverse momentum, JHEP 10 (2016) 066 [arXiv:1607.08817] [INSPIRE].

[30] K. Melnikov, L. Tancredi and C. Wever, Two-loop amplitudes for $qg \to Hq$ and $q\bar{q} \to Hg$ mediated by a nearly massless quark, Phys. Rev. D 95 (2017) 054012 [arXiv:1702.00426] [INSPIRE].

[31] J.M. Lindert, K. Melnikov, L. Tancredi and C. Wever, Top-bottom interference effects in Higgs plus jet production at the LHC, Phys. Rev. Lett. 118 (2017) 252002 [arXiv:1703.03886] [INSPIRE].

[32] J.M. Lindert, K. Kudashkin, K. Melnikov and C. Wever, Higgs bosons with large transverse momentum at the LHC, Phys. Lett. B 782 (2018) 210 [arXiv:1801.08226] [INSPIRE].

[33] T. Neumann, NLO Higgs+jet at Large Transverse Momentum Including Top Quark Mass Effects, J. Phys. Comm. 2 (2018) 095017 [arXiv:1802.02981] [INSPIRE].

[34] F. Caola, J.M. Lindert, K. Melnikov, P.F. Monni, L. Tancredi and C. Wever, Bottom-quark effects in Higgs production at intermediate transverse momentum, JHEP 09 (2018) 035 [arXiv:1804.07632] [INSPIRE].

[35] S.P. Jones, M. Kerner and G. Luisoni, Next-to-Leading-Order QCD Corrections to Higgs Boson Plus Jet Production with Full Top-Quark Mass Dependence, Phys. Rev. Lett. 120 (2018) 162001 [arXiv:1802.00349] [INSPIRE].

[36] R. Hamberg, W.L. van Neerven and T. Matsuura, A complete calculation of the order $\alpha_s^2$ correction to the Drell-Yan K factor, Nucl. Phys. B 359 (1991) 343 [Erratum ibid. B 644 (2002) 403] [INSPIRE].

[37] W.L. van Neerven and E.B. Zijlstra, The $O(\alpha_s^2)$ corrected Drell-Yan K factor in the DIS and MS scheme, Nucl. Phys. B 382 (1992) 11 [Erratum ibid. B 680 (2004) 513] [INSPIRE].

[38] C. Anastasiou, L.J. Dixon, K. Melnikov and F. Petriello, Dilepton rapidity distribution in the Drell-Yan process at NNLO in QCD, Phys. Rev. Lett. 91 (2003) 182002 [hep-ph/0306192] [INSPIRE].
[39] K. Melnikov and F. Petriello, The W boson production cross section at the LHC through $O(\alpha_s^2)$, Phys. Rev. Lett. 96 (2006) 231803 [hep-ph/0603182] [INSPIRE].

[40] K. Melnikov and F. Petriello, Electroweak gauge boson production at hadron colliders through $O(\alpha_s^2)$, Phys. Rev. D 74 (2006) 114017 [hep-ph/0609070] [INSPIRE].

[41] S. Catani, G. Ferrera and M. Grazzini, W Boson Production at Hadron Colliders: The Lepton Charge Asymmetry in NNLO QCD, JHEP 05 (2010) 006 [arXiv:1002.3115] [INSPIRE].

[42] S. Catani, L. Cieri, G. Ferrera, D. de Florian and M. Grazzini, Vector boson production at hadron colliders: a fully exclusive QCD calculation at NNLO, Phys. Rev. Lett. 103 (2009) 082001 [arXiv:0903.2120] [INSPIRE].

[43] R. Gavin, Y. Li, F. Petriello and S. Quackenbush, FEWZ 2.0: A code for hadronic Z production at next-to-next-to-leading order, Comput. Phys. Commun. 182 (2011) 2388 [arXiv:1011.3540] [INSPIRE].

[44] C. Anastasiou, L.J. Dixon, K. Melnikov and F. Petriello, High precision QCD at hadron colliders: Electroweak gauge boson rapidity distributions at NNLO, Phys. Rev. D 69 (2004) 094008 [hep-ph/0312266] [INSPIRE].

[45] A. Gehrmann-De Ridder, T. Gehrmann, E.W.N. Glover, A. Huss and T.A. Morgan, Precise QCD predictions for the production of a Z boson in association with a hadronic jet, Phys. Rev. Lett. 117 (2016) 022001 [arXiv:1507.02850] [INSPIRE].

[46] A. Gehrmann-De Ridder, T. Gehrmann, E.W.N. Glover, A. Huss and T.A. Morgan, The NNLO QCD corrections to Z boson production at large transverse momentum, JHEP 07 (2016) 133 [arXiv:1605.04295] [INSPIRE].

[47] A. Gehrmann-De Ridder, T. Gehrmann, E.W.N. Glover, A. Huss and T.A. Morgan, NNLO QCD corrections for Drell-Yan $p_T^Z$ and $\phi^*$ observables at the LHC, JHEP 11 (2016) 094 [Erratum ibid. 1810 (2018) 126] [arXiv:1610.01843] [INSPIRE].

[48] R. Gauld, A. Gehrmann-De Ridder, T. Gehrmann, E.W.N. Glover and A. Huss, Precise predictions for the angular coefficients in Z-boson production at the LHC, JHEP 11 (2017) 003 [arXiv:1708.00008] [INSPIRE].

[49] R. Boughezal et al., Z-boson production in association with a jet at next-to-next-to-leading order in perturbative QCD, Phys. Rev. Lett. 116 (2016) 152001 [arXiv:1512.01291] [INSPIRE].

[50] R. Boughezal, X. Liu and F. Petriello, Phenomenology of the Z-boson plus jet process at NNLO, Phys. Rev. D 94 (2016) 074015 [arXiv:1602.08140] [INSPIRE].

[51] R. Boughezal, C. Focke, X. Liu and F. Petriello, W-boson production in association with a jet at next-to-next-to-leading order in perturbative QCD, Phys. Rev. Lett. 115 (2015) 062002 [arXiv:1504.02131] [INSPIRE].

[52] R. Boughezal, X. Liu and F. Petriello, W-boson plus jet differential distributions at NNLO in QCD, Phys. Rev. D 94 (2016) 113009 [arXiv:1602.06965] [INSPIRE].

[53] A. Gehrmann-De Ridder, T. Gehrmann, E.W.N. Glover, A. Huss and D.M. Walker, Next-to-Next-to-Leading-Order QCD Corrections to the Transverse Momentum Distribution of Weak Gauge Bosons, Phys. Rev. Lett. 120 (2018) 122001 [arXiv:1712.07543] [INSPIRE].
[54] J.H. Kühn, A. Kulesza, S. Pozzorini and M. Schulze, *One-loop weak corrections to hadronic production of Z bosons at large transverse momenta*, Nucl. Phys. B 727 (2005) 368 [hep-ph/0507178] [INSPIRE].

[55] J.H. Kühn, A. Kulesza, S. Pozzorini and M. Schulze, *Electroweak corrections to large transverse momentum production of W bosons at the LHC*, Phys. Lett. B 651 (2007) 160 [hep-ph/0703283] [INSPIRE].

[56] A. Denner, S. Dittmaier, T. Kasprzik and A. Mück, *Electroweak corrections to W + jet hadroproduction including leptonic W-boson decays*, JHEP 08 (2009) 075 [arXiv:0906.1656] [INSPIRE].

[57] A. Denner, S. Dittmaier, T. Kasprzik and A. Mück, *Electroweak corrections to dilepton + jet production at hadron colliders*, JHEP 06 (2011) 069 [arXiv:1103.0914] [INSPIRE].

[58] G. Parisi and R. Petronzio, *Small Transverse Momentum Distributions in Hard Processes*, Nucl. Phys. B 154 (1979) 427 [INSPIRE].

[59] J.C. Collins, D.E. Soper and G.F. Sterman, *Transverse Momentum Distribution in Drell-Yan Pair and W and Z Boson Production*, Nucl. Phys. B 250 (1985) 199 [INSPIRE].

[60] G. Bozzi, S. Catani, D. de Florian and M. Grazzini, *Transverse-momentum resummation and the spectrum of the Higgs boson at the LHC*, Nucl. Phys. B 737 (2006) 73 [hep-ph/0508068] [INSPIRE].

[61] D. de Florian, G. Ferrera, M. Grazzini and D. Tommasini, *Higgs boson production at the LHC: transverse momentum resummation effects in the $H \to \gamma\gamma$, $H \to WW \to l\bar{l}l\nu$ and $H \to ZZ \to 4l$ decay modes*, JHEP 06 (2012) 132 [arXiv:1203.6321] [INSPIRE].

[62] T. Becher, M. Neubert and D. Wilhelm, *Higgs-Boson Production at Small Transverse Momentum*, JHEP 05 (2013) 110 [arXiv:1212.2621] [INSPIRE].

[63] C.W. Bauer, S. Fleming and M.E. Luke, *Summing Sudakov logarithms in $\vec{B} \to X_s\gamma$ in effective field theory*, Phys. Rev. D 63 (2000) 014006 [hep-ph/0005275] [INSPIRE].

[64] C.W. Bauer, S. Fleming, D. Pirjol and I.W. Stewart, *An Effective field theory for collinear and soft gluons: Heavy to light decays*, Phys. Rev. D 63 (2001) 114020 [hep-ph/0011336] [INSPIRE].

[65] C.W. Bauer, D. Pirjol and I.W. Stewart, *Soft collinear factorization in effective field theory*, Phys. Rev. D 65 (2002) 054022 [hep-ph/0109045] [INSPIRE].

[66] C.W. Bauer, S. Fleming, D. Pirjol, I.Z. Rothstein and I.W. Stewart, *Hard scattering factorization from effective field theory*, Phys. Rev. D 66 (2002) 014017 [hep-ph/0202088] [INSPIRE].

[67] T. Becher and M. Neubert, *Drell-Yan Production at Small $q_T$, Transverse Parton Distributions and the Collinear Anomaly*, Eur. Phys. J. C 71 (2011) 1665 [arXiv:1007.4005] [INSPIRE].

[68] M.G. Echevarria, A. Idilbi and I. Scimemi, *Factorization Theorem For Drell-Yan At Low $q_T$ And Transverse Momentum Distributions On-The-Light-Cone*, JHEP 07 (2012) 002 [arXiv:1111.4996] [INSPIRE].

[69] D. Neill, I.Z. Rothstein and V. Vaidya, *The Higgs Transverse Momentum Distribution at NNLL and its Theoretical Errors*, JHEP 12 (2015) 097 [arXiv:1503.00005] [INSPIRE].
[70] G. Bozzi, S. Catani, G. Ferrera, D. de Florian and M. Grazzini, Production of Drell-Yan lepton pairs in hadron collisions: Transverse-momentum resummation at next-to-next-to-leading logarithmic accuracy, Phys. Lett. B 696 (2011) 207 [arXiv:1007.2351] [INSPIRE].

[71] A. Banfi, M. Dasgupta, S. Marzani and L. Tomlinson, Predictions for Drell-Yan $\phi^*$ and $Q_T$ observables at the LHC, Phys. Lett. B 715 (2012) 152 [arXiv:1205.4760] [INSPIRE].

[72] H.-n. Li, Unification of the $k_T$ and threshold resummations, Phys. Lett. B 454 (1999) 328 [hep-ph/9812363] [INSPIRE].

[73] E. Laenen, G.F. Sterman and W. Vogelsang, Recoil and threshold corrections in short distance cross-sections, Phys. Rev. D 63 (2001) 114018 [hep-ph/0010080] [INSPIRE].

[74] A. Kulesza, G.F. Sterman and W. Vogelsang, Joint resummation for Higgs production, Phys. Rev. D 69 (2004) 014012 [hep-ph/0309264] [INSPIRE].

[75] S. Marzani, Combining $Q_T$ and small-x resummations, Phys. Rev. D 93 (2016) 054047 [arXiv:1511.06039] [INSPIRE].

[76] S. Forte and C. Muselli, High energy resummation of transverse momentum distributions: Higgs in gluon fusion, JHEP 03 (2016) 122 [arXiv:1511.05561] [INSPIRE].

[77] F. Caola, S. Forte, S. Marzani, C. Muselli and G. Vita, The Higgs transverse momentum spectrum with finite quark masses beyond leading order, JHEP 08 (2016) 150 [arXiv:1606.04100] [INSPIRE].

[78] G. Lustermans, W.J. Waalewijn and L. Zeune, Joint transverse momentum and threshold resummation beyond NLL, Phys. Lett. B 762 (2016) 447 [arXiv:1605.02740] [INSPIRE].

[79] S. Marzani and V. Theeuwes, Vector boson production in joint resummation, JHEP 02 (2017) 127 [arXiv:1612.01432] [INSPIRE].

[80] C. Muselli, S. Forte and G. Ridolfi, Combined threshold and transverse momentum resummation for inclusive observables, JHEP 03 (2017) 106 [arXiv:1701.01464] [INSPIRE].

[81] R.K. Ellis and S. Veseli, W and Z transverse momentum distributions: Resummation in $q_T$ space, Nucl. Phys. B 511 (1998) 649 [hep-ph/9706526] [INSPIRE].

[82] S. Frixione, P. Nason and G. Ridolfi, Problems in the resummation of soft gluon effects in the transverse momentum distributions of massive vector bosons in hadronic collisions, Nucl. Phys. B 542 (1999) 311 [hep-ph/9809367] [INSPIRE].

[83] A. Kulesza and W.J. Stirling, On the resummation of subleading logarithms in the transverse momentum distribution of vector bosons produced at hadron colliders, JHEP 01 (2000) 016 [hep-ph/9909271] [INSPIRE].

[84] P.F. Monni, E. Re and P. Torrielli, Higgs Transverse-Momentum Resummation in Direct Space, Phys. Rev. Lett. 116 (2016) 242001 [arXiv:1604.02191] [INSPIRE].

[85] W. Bizon, P.F. Monni, E. Re, L. Rottoli and P. Torrielli, Momentum-space resummation for transverse observables and the Higgs $p_\perp$ at $N^3LL+NNLO$, JHEP 02 (2018) 108 [arXiv:1705.09127] [INSPIRE].

[86] M.A. Ebert and F.J. Tackmann, Resummation of transverse momentum distributions in distribution space, JHEP 02 (2017) 110 [arXiv:1611.08610] [INSPIRE].
[87] D. Kang, C. Lee and V. Vaidya, *A fast and accurate method for perturbative resummation of transverse momentum-dependent observables*, JHEP 04 (2018) 149 [arXiv:1710.00078] [INSPIRE].

[88] S. Catani and M. Grazzini, *Higgs boson production at hadron colliders: hard-collinear coefficients at the NNLO*, Eur. Phys. J. C 72 (2012) 2013 [Erratum ibid. C 72 (2012) 2132] [arXiv:1106.4652] [INSPIRE].

[89] S. Catani, L. Cieri, D. de Florian, G. Ferrera and M. Grazzini, *Vector boson production at hadron colliders: hard-collinear coefficients at the NNLO*, Eur. Phys. J. C 72 (2012) 2195 [arXiv:1209.0158] [INSPIRE].

[90] T. Gehrmann, T. Luebbert and L.L. Yang, *Calculation of the transverse parton distribution functions at next-to-next-to-leading order*, JHEP 06 (2014) 155 [arXiv:1403.6451] [INSPIRE].

[91] M.G. Echevarria, I. Scimemi and A. Vladimirov, *Unpolarized Transverse Momentum Dependent Parton Distribution and Fragmentation Functions at next-to-next-to-leading order*, JHEP 09 (2016) 004 [arXiv:1604.07869] [INSPIRE].

[92] Y. Li and H.X. Zhu, *Bootstrapping rapidity anomalous dimensions for transverse-momentum resummation*, Phys. Rev. Lett. 118 (2017) 022004 [arXiv:1604.01404] [INSPIRE].

[93] A.A. Vladimirov, *Correspondence between soft and rapidity anomalous dimensions*, Phys. Rev. Lett. 118 (2017) 062001 [arXiv:1610.05791] [INSPIRE].

[94] S. Moch, B. Ruijl, T. Ueda, J.A.M. Vermaseren and A. Vogt, *Four-loop non-singlet splitting functions in the planar limit and beyond*, JHEP 10 (2017) 041 [arXiv:1707.08315] [INSPIRE].

[95] S. Moch, B. Ruijl, T. Ueda, J.A.M. Vermaseren and A. Vogt, *On quartic colour factors in splitting functions and the gluon cusp anomalous dimension*, Phys. Lett. B 782 (2018) 627 [arXiv:1805.09638] [INSPIRE].

[96] X. Chen et al., *Precise QCD Description of the Higgs Boson Transverse Momentum Spectrum*, Phys. Lett. B 788 (2019) 425 [arXiv:1805.00736] [INSPIRE].

[97] F. Wilczek, *Decays of Heavy Vector Mesons Into Higgs Particles*, Phys. Rev. Lett. 39 (1977) 1304 [INSPIRE].

[98] M.A. Shifman, A.I. Vainshtein and V.I. Zakharov, *Remarks on Higgs Boson Interactions with Nucleons*, Phys. Lett. B 78 (1978) 443 [INSPIRE].

[99] T. Inami, T. Kubota and Y. Okada, *Effective Gauge Theory and the Effect of Heavy Quarks in Higgs Boson Decays*, Z. Phys. C 18 (1983) 69 [INSPIRE].

[100] K.G. Chetyrkin, B.A. Kniehl and M. Steinhauser, *Decoupling relations to $O(\alpha_S^3)$ and their connection to low-energy theorems*, Nucl. Phys. B 510 (1998) 61 [hep-ph/9708255] [INSPIRE].

[101] ATLAS collaboration, *Measurement of the transverse momentum and $\phi_\eta^*$ distributions of Drell-Yan lepton pairs in proton-proton collisions at $\sqrt{s} = 8$ TeV with the ATLAS detector*, Eur. Phys. J. C 76 (2016) 291 [arXiv:1512.02192] [INSPIRE].

[102] V. Del Duca, A. Frizzo and F. Maltoni, *Higgs boson production in association with three jets*, JHEP 05 (2004) 064 [hep-ph/0404013] [INSPIRE].
[103] L.J. Dixon, E.W.N. Glover and V.V. Khoze, MHV rules for Higgs plus multi-gluon amplitudes, JHEP 12 (2004) 015 [hep-th/0411092] [INSPIRE].

[104] S.D. Badger, E.W.N. Glover and V.V. Khoze, MHV rules for Higgs plus multi-parton amplitudes, JHEP 03 (2005) 023 [hep-th/0412275] [INSPIRE].

[105] J.M. Campbell and R.K. Ellis, Next-to-leading order corrections to $W + 2$ jet and $Z + 2$ jet production at hadron colliders, Phys. Rev. D 65 (2002) 113007 [hep-ph/0202176] [INSPIRE].

[106] J.M. Campbell, R.K. Ellis and D.L. Rainwater, Next-to-leading order QCD predictions for $W + 2$ jet and $Z + 2$ jet production at the CERN LHC, Phys. Rev. D 68 (2003) 094021 [hep-ph/0308195] [INSPIRE].

[107] K. Hagiwara and D. Zeppenfeld, Amplitudes for Multiparton Processes Involving a Current at $e^{+}e^{-}$, $e^{\pm}p$ and Hadron Colliders, Nucl. Phys. B 313 (1989) 560 [INSPIRE].

[108] F.A. Berends, W.T. Giele and H. Kuijf, Exact Expressions for Processes Involving a Vector Boson and Up to Five Partons, Nucl. Phys. B 321 (1989) 39 [INSPIRE].

[109] N.K. Falck, D. Graudenz and G. Kramer, Cross-section for Five Jet Production in $e^{+}e^{-}$ Annihilation, Nucl. Phys. B 328 (1989) 317 [INSPIRE].

[110] L.J. Dixon and Y. Sofianatos, Analytic one-loop amplitudes for a Higgs boson plus four partons, JHEP 08 (2009) 058 [arXiv:0906.0008] [INSPIRE].

[111] S. Badger, E.W. Nigel Glover, P. Mastrolia and C. Williams, One-loop Higgs plus four gluon amplitudes: Full analytic results, JHEP 01 (2010) 036 [arXiv:0909.4475] [INSPIRE].

[112] S. Badger, J.M. Campbell, R.K. Ellis and C. Williams, Analytic results for the one-loop NMHV Hqgg amplitude, JHEP 12 (2009) 035 [arXiv:0910.4481] [INSPIRE].

[113] E.W.N. Glover and D.J. Miller, The One loop QCD corrections for $\gamma^{*} \to Q\bar{Q}q\bar{q}$, Phys. Lett. B 396 (1997) 257 [hep-ph/9609474] [INSPIRE].

[114] Z. Bern, L.J. Dixon, D.A. Kosower and S. Weinzierl, One loop amplitudes for $e^{+}e^{-} \to \bar{q}q\bar{Q}Q$, Nucl. Phys. B 489 (1997) 3 [hep-ph/9610370] [INSPIRE].

[115] J.M. Campbell, E.W.N. Glover and D.J. Miller, The One loop QCD corrections for $\gamma^{*} \to q\bar{q}gg$, Phys. Lett. B 409 (1997) 503 [hep-ph/9706297] [INSPIRE].

[116] Z. Bern, L.J. Dixon and D.A. Kosower, One loop amplitudes for $e^{+}e^{-}$ to four partons, Nucl. Phys. B 513 (1998) 3 [hep-ph/9708239] [INSPIRE].

[117] T. Gehrmann, M. Jaquier, E.W.N. Glover and A. Koukoutsakis, Two-Loop QCD Corrections to the Helicity Amplitudes for $H \to 3$ partons, JHEP 02 (2012) 056 [arXiv:1112.3554] [INSPIRE].

[118] S. Moch, P. Uwer and S. Weinzierl, Two loop amplitudes with nested sums: Fermionic contributions to $e^{+}e^{-} \to q\bar{q}g$, Phys. Rev. D 66 (2002) 114001 [hep-ph/0207043] [INSPIRE].

[119] L.W. Garland, T. Gehrmann, E.W.N. Glover, A. Koukoutsakis and E. Remiddi, The Two loop QCD matrix element for $e^{+}e^{-} \to 3$ jets, Nucl. Phys. B 627 (2002) 107 [hep-ph/0112081] [INSPIRE].

[120] L.W. Garland, T. Gehrmann, E.W.N. Glover, A. Koukoutsakis and E. Remiddi, Two loop QCD helicity amplitudes for $e^{+}e^{-} \to 3$ jets, Nucl. Phys. B 642 (2002) 227 [hep-ph/0206067] [INSPIRE].

[121] T. Gehrmann and L. Tancredi, Two-loop QCD helicity amplitudes for $q\bar{q} \to W^{\pm}\gamma$ and $q\bar{q} \to Z^{0}\gamma$, JHEP 02 (2012) 004 [arXiv:1112.1531] [INSPIRE].
[122] A. Gehrmann-De Ridder, T. Gehrmann and E.W.N. Glover, quark-gluon antenna functions from neutralino decay, Phys. Lett. B 612 (2005) 36 [hep-ph/0501291] [INSPIRE].

[123] A. Daleo, T. Gehrmann and D. Maître, Antenna subtraction with hadronic initial states, JHEP 04 (2007) 016 [hep-ph/0612257] [INSPIRE].

[124] J. Currie, E.W.N. Glover and S. Wells, Infrared Structure at NNLO Using Antenna Subtraction, JHEP 04 (2013) 066 [arXiv:1301.4693] [INSPIRE].

[125] X. Chen, T. Gehrmann, E.W.N. Glover and M. Jaquier, Precise QCD predictions for the production of Higgs + jet final states, Phys. Lett. B 740 (2015) 147 [arXiv:1408.5325] [INSPIRE].

[126] A. Gehrmann-De Ridder, T. Gehrmann and E.W.N. Glover, Antenna subtraction at NNLO, JHEP 09 (2005) 056 [hep-ph/0505111] [INSPIRE].

[127] A. Gehrmann-De Ridder, T. Gehrmann and E.W.N. Glover, Gluon-gluon antenna functions from Higgs boson decay, Phys. Lett. B 612 (2005) 49 [hep-ph/0502110] [INSPIRE].

[128] A. Daleo, A. Gehrmann-De Ridder, T. Gehrmann and G. Luisoni, Antenna subtraction at NNLO with hadronic initial states: initial-final configurations, JHEP 01 (2010) 118 [arXiv:0912.0374] [INSPIRE].

[129] R. Boughezal, A. Gehrmann-De Ridder and M. Ritzmann, Antenna subtraction at NNLO with hadronic initial states: double real radiation for initial-initial configurations with two quark flavours, JHEP 02 (2011) 098 [arXiv:1011.6631] [INSPIRE].

[130] T. Gehrmann and P.F. Monni, Antenna subtraction at NNLO with hadronic initial states: real-virtual initial-initial configurations, JHEP 12 (2011) 049 [arXiv:1107.4037] [INSPIRE].

[131] A. Gehrmann-De Ridder, T. Gehrmann and M. Ritzmann, Antenna subtraction at NNLO with hadronic initial states: double real initial-initial configurations, JHEP 10 (2012) 047 [arXiv:1207.5779] [INSPIRE].

[132] L.J. Dixon, L. Magnea and G.F. Sterman, Universal structure of subleading infrared poles in gauge theory amplitudes, JHEP 08 (2008) 022 [arXiv:0805.3515] [INSPIRE].

[133] A. Banfi, G.P. Salam and G. Zanderighi, Principles of general final-state resummation and automated implementation, JHEP 03 (2005) 073 [hep-ph/0407286] [INSPIRE].

[134] A. Banfi, H. McAslan, P.F. Monni and G. Zanderighi, A general method for the resummation of event-shape distributions in $e^{+}e^{-}$ annihilation, JHEP 05 (2015) 102 [arXiv:1412.2126] [INSPIRE].

[135] S. Moch, J.A.M. Vermaseren and A. Vogt, Higher-order corrections in threshold resummation, Nucl. Phys. B 726 (2005) 317 [hep-ph/0506288] [INSPIRE].

[136] S. Moch, J.A.M. Vermaseren and A. Vogt, The three loop splitting functions in QCD: The nonsinglet case, Nucl. Phys. B 688 (2004) 101 [hep-ph/0403192] [INSPIRE].

[137] A. Vogt, S. Moch and J.A.M. Vermaseren, The three-loop splitting functions in QCD: The singlet case, Nucl. Phys. B 691 (2004) 129 [hep-ph/0404111] [INSPIRE].

[138] A. Banfi, G.P. Salam and G. Zanderighi, NLL+NNLO predictions for jet-veto efficiencies in Higgs-boson and Drell-Yan production, JHEP 06 (2012) 159 [arXiv:1203.5773] [INSPIRE].

[139] A. Banfi, P.F. Monni, G.P. Salam and G. Zanderighi, Higgs and Z-boson production with a jet veto, Phys. Rev. Lett. 109 (2012) 202001 [arXiv:1206.4998] [INSPIRE].
[140] A. Banfi et al., Jet-vetoed Higgs cross section in gluon fusion at N³LO+NNLL with small-R resummation, JHEP 04 (2016) 049 [arXiv:1511.02886] [INSPIRE].

[141] NNPDF collaboration, R.D. Ball et al., Parton distributions for the LHC Run II, JHEP 04 (2015) 040 [arXiv:1410.8849] [INSPIRE].

[142] J. Butterworth et al., PDF4LHC recommendations for LHC Run II, J. Phys. G 43 (2016) 023001 [arXiv:1510.03865] [INSPIRE].

[143] S. Dulat et al., New parton distribution functions from a global analysis of quantum chromodynamics, Phys. Rev. D 93 (2016) 033006 [arXiv:1506.07443] [INSPIRE].

[144] L.A. Harland-Lang, A.D. Martin, P. Motylinski and R.S. Thorne, Parton distributions in the LHC era: MMHT 2014 PDFs, Eur. Phys. J. C 75 (2015) 204 [arXiv:1412.3989] [INSPIRE].

[145] S. Carrazza, J.I. Latorre, J. Rojo and G. Watt, A compression algorithm for the combination of PDF sets, Eur. Phys. J. C 75 (2015) 474 [arXiv:1504.06469] [INSPIRE].

[146] G. Watt and R.S. Thorne, Study of Monte Carlo approach to experimental uncertainty propagation with MSTW 2008 PDFs, JHEP 08 (2012) 052 [arXiv:1205.4024] [INSPIRE].

[147] P. Pietrulewicz, D. Samitz, A. Spiering and F.J. Tackmann, Factorization and Resummation for Massive Quark Effects in Exclusive Drell-Yan, JHEP 08 (2017) 114 [arXiv:1703.09702] [INSPIRE].

[148] E. Bagnaschi, F. Maltoni, A. Vicini and M. Zaro, Lepton-pair production in association with a $b\bar{b}$ pair and the determination of the W boson mass, JHEP 07 (2018) 101 [arXiv:1803.04336] [INSPIRE].

[149] D. de Florian and M. Grazzini, The Structure of large logarithmic corrections at small transverse momentum in hadronic collisions, Nucl. Phys. B 616 (2001) 247 [hep-ph/0108273] [INSPIRE].

[150] C.T.H. Davies and W.J. Stirling, Nonleading Corrections to the Drell-Yan Cross-Section at Small Transverse Momentum, Nucl. Phys. B 244 (1984) 337 [INSPIRE].