Fiducial distributions in Higgs and Drell-Yan production at N$^3$LL+NNLO

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.


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][2][3][4][5][6][7][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 transversemomentum-dependent PDFs [10][11][12]). An example is given by the differential distributions in Zand 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 3 LL), for the transversemomentum spectrum of the colour singlet, as well as for the angular variable φ * η [17]. In the following, for simplicity, we will collectively denote p t /M or φ * η 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-tonext-to-leading-order (N 3 LO) accuracy in QCD [18][19][20][21][22][23][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][27][28][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][31][32][33][34][35].
Although fixed-order results are crucial to obtain reliable theoretical predictions away from the soft and collinear regions of the phase space (v ∼ 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 α n s ln k (1/v)/v, with k ≤ 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 Σ as ln Σ(v) ≡ ln the small-transverse-momentum region was also studied in detail in Refs. [72][73][74][75][76][77][78][79][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][82][83], but remained unsolved until recently. Due to the vectorial nature of p t (analogous considerations apply to φ * η ), 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 3 LL resummation of p t (and φ * η ) spectra in color-singlet production have been computed in [88][89][90][91][92], with the exception of the four-loop cusp anomalous dimension, which is currently known only for leading colour contributions [93]. 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 3 LL matched to NNLO in Refs. [85,94]. In this manuscript, employing the direct-space resummation at N 3 LL 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 → γγ channel. We also consider Drell-Yan pair production and compute N 3 LL+NNLO predictions for the transverse momentum of the lepton pair and for the φ * η 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 φ * η 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 3 LL.

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 φ * η 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 H t 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 → ∞), giving rise to an effective operator that directly couples the Higgs field to the gluon field-strength tensor via [95][96][97] The Wilson coefficient λ is known to three-loop accuracy [98] and its renormalisation-scale dependence was studied in [29]. We consider the p H t 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 [99], which are summarised in Section 6. We consider both the p Z t spectrum as well as the φ * η distribution, which are further studied multi-differentially for different invariant-mass (M ) or rapidity (Y ) bins.
The differential distributions in v = p t /M, φ * η 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 + jet process with X = H, Z, where the jet requirement is replaced by a restriction on v to be non-vanishing: v ≥ v cut > 0. The state-of-the-art fixed-order QCD predictions for this class of processes is at NNLO [26][27][28][29][45][46][47][48][49][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: doublereal corrections for H [100][101][102] and DY [103][104][105][106][107]), with one extra parton and one extra loop (RV: real-virtual corrections for H [108][109][110] and DY [103,104,[111][112][113][114]) and with no extra parton but two extra loops (VV: double-virtual corrections for H [115] and DY [116][117][118][119]). 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 a 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 [120][121][122] 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 → H +jet [29,123] and pp → Z +jet [45][46][47][48] at NNLO, and it takes the schematic form: The antenna subtraction terms, dσ S,T,U NNLO , for both Higgs and Drell-Yan related processes are constructed from antenna functions [124][125][126][127][128][129] to cancel infrared singularities between the contributions of different parton multiplicities. The integrals are performed over the phase space Φ 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σ NNLO X+jet /dO. For large values of v (v ∼ 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 (∼ 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 → 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 φ * η = 0.004 [47]. The accuracy of the results obtained with the NNLOjet code for small v has been systematically validated in Ref. [94] by comparing fixed-order predictions of the Higgs boson transverse momentum distribution dσ NNLO /dp H t against the expansion of the N 3 LL resummation (obtained in the framework of soft-collinear effective field theory, SCET) to the respective order in the small p H t region. This validation was performed for individual initial-state partonic channels down to p H t = 0.7 GeV.
As v → 0, the final-state phase space Φ X+1,2,3 is reduced to the phase space of colour singlet production Φ X . The RR, RV, and VV contributions contain infrared divergences with one extra unresolved parton that cannot be cancelled by the subtraction terms dσ S,T,U NNLO . 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.

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 for an observable v ( ) = V (Φ B , k 1 , . . . , k n ), being either p t /M or φ * η , in the presence of n real emissions with momenta k 1 , ..., k n . Σ(v) can be expressed as where M is the matrix element for n real emissions and V(Φ B ) denotes the resummed form factor that encodes the purely virtual corrections [130]. The phase spaces of the i-th emission k i and that of the Born configuration 1 are denoted by [dk i ] and dΦ B , respectively. The recursive infrared and collinear (rIRC) safety [131] of the observable allows one to establish a well defined logarithmic counting in the squared amplitude [131,132], 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 V(Φ 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 (Φ B , k 1 , . . . , k n ), thus it exponentiates and cancels the divergences contained in V(Φ 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 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 Θ 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 → 0 can be taken safely.
For observables which depend on the total transverse momentum of QCD radiation, such as p t or φ * η , it is particularly convenient to set the resolution scale to a small fraction δ > 0 of the transverse momentum of the block with largest k t , hereby denoted by δ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 φ * η . Including terms up to N 3 LL, the cumulative cross section in momentum space can be recast in the following form [85] where ζ si ≡ k tsi /k t1 and we introduced the notation dZ[{R , k i }] to denote an ensemble that describes the emission of n identical independent blocks [85]. The average of a function G({p}, {k i }) over the measure dZ is defined as ( (3.4) The ln 1/δ 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 δ-independent.
To obtain Eq. (3.3) we used the fact that, for resolved radiation, ζ 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 k t1 do not contribute to the value of the observable and are therefore fully cancelled by the term exp{−R (k t1 ) ln(1/δ)} of Eq. (3.4). Although not strictly necessary, this expansion allows for a more efficient numerical implementation. The expansion gives rise to the termsR (n) which denote the derivatives of the radiator as whereR takes the form and α s = α s (µ R ). We report the functions g i in Appendix B, and we refer to Ref. [85] for further details.
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 Q. Such modified logarithmsL 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 → Q onto k t1 → ∞ in all terms of Eq. (3.3), with the exception of the observable's measurement function. A convenient choice of such a mapping is 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 J , defined as This corresponds to the Jacobian for the transformation (3.7), and ensures the absence of fractional (although power suppressed) α 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 Θ 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 factorsL contain the parton luminosities up to N 3 LL, multiplied by the Born-level squared, and virtual amplitudes. They are defined as Y is the rapidity of the colour singlet in the centre-of-mass frame of the collision at the Born-level, |M B | 2 cc is the Born-level squared matrix element, and x Q = Q/M . The above luminosities contain the NLO and NNLO coefficient functionsC ci for Higgs and Drell-Yan production [88][89][90], as well as the hard virtual correctionsH (n) . A precise definition is given is Section 4 of Ref. [85], and the relevant formulae are also reported in Appendix B.
Finally, we define the convolution of a regularised splitting functionP [133,134] with the coefficientL NLL aŝ The termP (0) ⊗P (0) ⊗L 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 . (3.14)

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 Σ, we define the analogue of Eq. (3.1) for the fixed-order prediction as tot is the total cross section for the considered processes and dΣ NNLO /dv denotes the NNLO differential distribution.
For inclusive Higgs production, the transverse-momentum distribution at NNLO was obtained in Refs. [26][27][28][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 → γγ with fiducial cuts, we approximate the N 3 LO correction to σ N 3 LO tot 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 where Σ EXP denotes the expansion of the resummation formula Σ N 3 LL to N 3 LO. The second scheme we consider belongs to the class of multiplicative schemes similar to those defined in Refs. [135][136][137], and it is schematically defined as where the quantity in square brackets is expanded to N 3 LO. The two schemes (4.2), (4.3) are equivalent at the perturbative order we are working at, and differ by N 4 LO and N 4 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 3 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 3 LO constant terms, of formal N 4 LL accuracy, are automatically extracted from the fixed order in the matching procedure, whenever the N 3 LO total cross section is known. We recall that Eq. (3.3) resums all towers of ln(1/v) up to N 3 LL, defined at the level of the logarithm of Σ (1.1). At this order, one predicts correctly all logarithmic terms up to, and including, α n s ln 2n−5 (1/v) in the expanded formula for Σ, while terms of order α n s ln 2n−6 (1/v) would be modified by including N 4 LL corrections. The inclusion of constant terms of order O(α 3 s ) relative to Born level in the resummed formula, of formal N 4 LL accuracy, extends the prediction to all terms of order α n s ln 2n−6 (1/v) in the expanded formula for Σ. Indeed these terms, which contain the N 3 LO collinear coefficient functions and three-loop virtual corrections, would multiply the Sudakov e −R(kt1) in the resummed formula (3.3) starting at N 4 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 4 LL terms α n s ln 2n−6 (1/v) in the expanded formula for the matched cumulative cross section Σ MAT . 3 We stress that this, as pointed out above, requires the knowledge of the N 3 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 3 LO cross sections are approximated, as described at the beginning of this section, and hence the tower of N 4 LL terms α n s 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 limitL → 0, Σ N 3 LL tends to the integral ofL N 3 LL (µ F ) (defined in Eq. (3.11)) over Φ B , evaluated atL = 0. Therefore, the fixed-order result Σ N 3 LO at large v receives a spurious correction of relative order α 4 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 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 L → 0 value. This is simply given by the integral over the Born phase space Φ B of theL → 0 limit ofL N 3 LL (that we report in Eq. (A.5)) where the integration over Φ B is performed by taking into account the phase-space cuts of the experimental analysis.
We thus obtain  and the whole squared bracket in Eq. (4.7) is expanded to N 3 LO. This ensures that, in the v 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 Σ N 3 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.

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 δ to zero, since the cancellation of IRC divergences is manifest. In Figure 1 we show the comparison between the expansion of the N 3 LL resummed cross section and the fixed order for the differential distribution of p H t 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 3 LO cumulant, and similarly for lower orders.
In Figure 1 we see that below p H t ∼ 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 H t ∼ 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 H t . In particular, we consider the differential in order to highlight potential logarithmic differences in the p H t → 0 region. A similar validation of the NNLO p H t distribution has been performed in Ref. [94]. 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(α 3 s ) expansion of the NNLL resummation. As one expects, at small p H t 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 3 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 3 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 H t 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 3 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 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.

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 3 LL resummed calculation. We use the distribution of the φ * η observable for the validation which, also from the experimental viewpoint, is particularly suitable to probe the deep infrared regime of the radiation phase space. This is done without loss of generality, as this distribution shares the same infrared behaviour with the p Z t distribution [138]. In particular, we focus on the differential distribution , (4. 10) in order to highlight potential logarithmic differences in the φ * η → 0 region. To perform the validation we consider 8 TeV pp collisions with NNPDF3.0 parton densities [139], and we work within the following fiducial cuts (corresponding to that used by the ATLAS collaboration [99]) where the large invariant mass is chosen in order to enhance the impact of logarithmic corrections. This slightly facilitates the comparison to fixed-order in the limit φ * η → 0. We set the scales as µ R = µ F = M T = M 2 + p 2 t , and x Q = Q/M = 1. The results 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, and between the NNLO distribution and the expansion of the N 3 LL resummation to third order. In both cases one expects the differences to approach zero at small φ * η , which is well confirmed by the plot. In general, the asymptotic region for DY production is shifted towards smaller p t (and hence φ * η ) values with respect to the inclusive Higgs production case that we studied in Section 4.1. This difference is partly driven by the smaller colour factor associated with the initial-state quarks, compared to  gluons in the Higgs case. Moreover, a further complication is given by the use of the 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, this makes it necessary to go to very small values of φ * η in order to see a convergence of the fixed-order to the expansion of the resummation, as shown in left panel of Figure 4. In order to single out the contribution of the NNLO correction, in the right panel of Figure 4 we show 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 φ * η values. 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.

Results for Higgs production in HEFT
In this section we present our predictions for the p H t spectrum both in inclusive pp → H production, and in the pp → H → γγ 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 [139][140][141][142][143][144]. The value of the parameter p appearing in the definition of the modified logarithmsL 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.
We set the central renormalisation and factorisation scales as µ R = µ 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 µ R , µ F by a factor of two in either direction, while keeping 1/2 < µ R /µ F < 2 and x Q = 1/2; Moreover, for central µ R and µ 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.

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 5 we compare the differential distributions at N 3 LL+NLO and NNLL+NLO in the small-p H t 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 H t . 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 H t 5 GeV. The right plot of Figure 5 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 H t . 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 ∼ 1/p H t at small p H t . The fact that the divergence is not visible in the distribution reported in the upper panel of Figure 5 is entirely due to the nature of the multiplicative scheme, which ensures that the distribution follows the resummation scaling at small p H t , 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 H t 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 H t , 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. [94].
Next, we consider the comparison between the matched prediction and the fixed-order one. Figure 6 shows this comparison for two different central scales. The left plot is obtained with central µ F = µ R = m H /2, while the right plot is obtained with µ F = µ R = m H . The rest of the setup is kept as described above. We observe that at µ F = µ 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 6) 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 H t 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.
To conclude, Figure 7 reports the comparison between our best prediction (N 3 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.

Matched predictions for fiducial H → γγ
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 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 → H → γγ and, in particular, we focus on the transverse momentum of the γγ system in the presence of fiducial cuts.
The fiducial volume is defined by the set of cuts detailed below where p γ1 t , p γ2 t are the transverse momenta of the two photons, η γi are their pseudo-rapidities in the hadronic centre-of-mass frame, and Y γγ 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 3 LL+NNLO accuracy of the differential distributions. 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 × 10 −3 .
In Figure 8 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 3 LL+NLO vs. NLO in the left panel, and N 3 LL+NNLO vs. NNLO in the right one.
By comparing the two panels of Figure 8 we notice a substantial reduction in the theoretical uncertainty in the medium-high-p γγ t region, driven by the increase in perturbative accuracy of the fixed-order computation; at very low p γγ t , 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 6. We 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 µ R = µ F = m H /2, which suffers from large accidental cancellations.

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 φ * η . We consider 8 TeV proton-proton collisions, and compare the resulting calculation for the differential spectra with ATLAS data from Ref. [99]. The fiducial phase-space volume is defined as follows: where p ± t are the transverse momenta of the two leptons, η ± are their pseudo-rapidities, while Y and M 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 µ R = µ F = M T = M 2 + p 2 t , while the central resummation scale is chosen to be x Q = Q/M = 1/2. The theoretical uncertainty is estimated through the same set of variations as for Higgs boson production.
The results for p Z t and φ * η are shown in the following two subsections. All plots have the same pattern: the main panels display the comparison of normalised differential distributions at NNLO (green), NNLL+NLO (blue), and N 3 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 Z t distributions In Figure 9 we display the normalised p Z t distributions in which, in addition to the fiducial cuts reported above, we consider three different lepton-pair invariant-mass windows: A comparison of the most accurate matched prediction with the fixed-order one shows that the N 3 LL+NNLO prediction starts differing significantly from the NNLO for p Z t < 10-15 GeV, while for p Z t > 20 GeV the NNLO is sufficient to provide a reliable description. Comparing matched predictions with different formal accuracy, we note that the N 3 LL+NNLO curve has a significantly reduced theoretical systematics with respect to that for NNLL+NLO, in the whole p Z t range and for all considered invariant-mass windows. The perturbative error is reduced by more than a factor of two at very low p Z t , 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 Z t distributions is also significantly distorted by the inclusion of higher orders: the spectrum is harder than the NNLL+NLO result for p Z t 10 GeV, and the peak is lower, with the N 3 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 3 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 Z t range. In Figure 10 we focus our analysis on the central lepton-pair invariant-mass window defined in Eq. (6.2) and show predictions for the normalised p Z t distribution in six different lepton-pair rapidity slices: The comments relevant to Figure 9 by far and large apply in this case as well, with our best prediction at N 3 LL+NNLO affected by an uncertainty that is of order 3-5% in the whole p Z t 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 Z t range, and the pure NNLO, in the p Z t 15 GeV region.
The pattern of comparisons among theoretical predictions is qualitatively similar to what discussed for the p Z t distribution. Resummation effects at N 3 LL+NNLO start being important with respect to the pure NNLO in the region φ * η 0.2; the shape of the N 3 LL+NNLO distribution is significantly distorted with respect to the NNLL+NLO one in a similar fashion as for the p Z t 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 Z t case, however, for φ * η we note that the N 3 LL+NNLO prediction describes data appropriately only in the central-and high-invariant-mass windows. In the lowinvariant-mass one, the prediction undershoots data in the medium-hard region, by up to 10%. This tension was already observed in the fixed-order NNLO comparison [47]. 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 φ * η , especially at low φ * η . In Figure 12 we show the results for the φ * η 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 Ratio to data Figure 9. Comparison of the normalised transverse momentum distribution for Drell-Yan pair production at NNLO (green), NNLL+NLO (blue) and N 3 LL+NNLO (red) at √ s = 8 TeV integrated over the full lepton-pair rapidity range (0 < |Y | < 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. availability of experimental measurements, in Figures 13 and 14 we also provide predictions sliced in Y 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).
The prediction subdivided in rapidity slices largely shares the same features as that integrated over rapidity, which has been detailed in Figure 11. 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 in most cases but for the low-φ * η region at large rapidity where the uncertainty is in the 5-10% range. 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. Ratio to data Figure 10. Comparison of the normalised transverse momentum distribution for Drell-Yan pair production at NNLO (green), NNLL+NLO (blue) and N 3 LL+NNLO (red) at √ s = 8 TeV in the central lepton-pair invariant-mass window (66 GeV < M < 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. Ratio to data Figure 11. Comparison of the normalised φ * η distribution for Drell-Yan pair production at NNLO (green), NNLL+NLO (blue) and N 3 LL+NNLO (red) at √ s = 8 TeV integrated over the full lepton-pair rapidity range (0 < |Y | < 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. Ratio to data Figure 12. Comparison of the normalised φ * η distribution for Drell-Yan pair production at NNLO (green), NNLL+NLO (blue) and N 3 LL+NNLO (red) at √ s = 8 TeV in the central lepton-pair invariant-mass window (66 GeV < M < 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. Ratio to data Figure 13. Comparison of the normalised φ * η distribution for Drell-Yan pair production at NNLO (green), NNLL+NLO (blue) and N 3 LL+NNLO (red) at √ s = 8 TeV in the low lepton-pair invariant-mass window (46 GeV < M < 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. Ratio to data Figure 14. Comparison of the normalised φ * η distribution for Drell-Yan pair production at NNLO (green), NNLL+NLO (blue) and N 3 LL+NNLO (red) at √ s = 8 TeV in the high lepton-pair invariant-mass window (116 GeV < M < 150 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.

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 → γγ 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 H t ∼ 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 H t . 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 [99] both for the normalised p Z t distributions and for the normalised φ * η 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 Z t < 10-15 GeV, while for p Z t > 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 Z t range and for all invariant-mass windows considered in our study. For the φ * η distribution, resummation effects start being important with respect to pure NNLO in the region φ * η 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 φ * η 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 Z t case, for φ * η 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 ∼ 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 H t 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 Z t and φ * η 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 [145,146], and more precise studies are necessary in order to assess their impact precisely. Recent analyses [146] 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 3 LL+NNLO prediction, a careful assessment will be necessary to improve further on the results presented in this work.

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 Moreover, we denote the perturbative expansion of the resummed cross section Σ N k LL as With this notation, the additive scheme of Eq. (4.2) becomes (for simplicity we drop the explicit dependence on v in the following) where the three terms in curly brackets denote the NLO, NNLO and N 3 LO contributions to the matching, respectively. For the multiplicative scheme we need to introduce the asymptotic expansion Σ N k LL asym. , defined in Eq. (4.6) (the definition for k = 3 is analogous with obvious replacements) in terms of theL → 0 limit of the coefficientsL N k LL of Eqs. (3.9), (3.10), (3.11), which read In the following formula the perturbative expansion of Σ N k LL asym. is denoted as follows With this notation the matching formula (4.7) reads asym.
where, as above, we grouped the terms entering at NLO, NNLO, and N 3 LO within curly brackets.

B Formulae for N 3 LL resummation
In this section we report the expressions for quantities needed for N 3 LL resummation of transverse observables, that we have used throughout this article.
We also provide expressions for the functions g i (λ) entering in the N 3 LL Sudakov radiator Eq. (3.6) and its derivative. We define (B.5) We have: (B.9) For Higgs boson production in gluon fusion, the coefficients A (i) and B (i) which enter the formulae above are (1) The expressions for the coefficients A (i) and B (i) are extracted from Refs. [62,91,92,147] for Higgs boson production and Refs. [67,91,92,148] for DY production. The N 3 LL anomalous dimension A (4) both for the Drell-Yan process and Higgs boson production is at present incomplete since the four-loop cusp anomalous dimension is unknown. In the above expressions, this contribution is set to zero. The new terms ∆A (4) = −64π 3 β 3 0 ζ 3 , ∆B (3) = −32π 2 β 2 0 ζ 3 , ∆H (2) = 16 3 πβ 0 ζ 3 , (B.12) are a feature due to performing the resummation in momentum space, and do not appear in the anomalous dimensions in b space (see Ref. [85] for details). The term ∆H (2) will appear in theH functions defined below. We also present the expansion of hard-virtual coefficient function H in powers of the strong coupling where µ is the same scale that enters parton densities. The first-order expansion has been known for a long time and reads The second-order collinear coefficient functions C (2) ab (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], 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 (2) gq and C (2) gg coefficients normalised as in Eq. (B.18) are extracted from Eqs. (30) and (32) of Ref. [88], respectively, where we use the hard coefficients of Eqs. (B.14) without the new term ∆H (2) in the H (2) g (M ) coefficient. 5 The coefficient G (1) is taken from Eq. (13) of Ref. [88]. Similarly, for quark-initiated processes, we extract C (2) qg and C (2) qq from Eqs. (32) and (34) of Ref. [89], respectively, where we use the hard coefficients from Eqs. (B.15) without the new term ∆H (2) in the H  [3] ATLAS collaboration, M. Aaboud et al., Measurement of differential cross sections and W + /W − cross-section ratios for W boson production in association with jets at √ s = 8 TeV with the ATLAS detector, 1711.03296.