Next-to-leading power threshold corrections for finite order and resummed colour-singlet cross sections

We study next-to-leading-power (NLP) threshold corrections in colour-singlet production processes, with particular emphasis on Drell-Yan (DY) and single-Higgs production. We assess the quality of the partonic and hadronic threshold expansions for each process up to NNLO. We determine numerically the NLP leading-logarithmic (LL) resummed contribution in addition to the leading-power next-to-next-to-leading logarithmic (LP NNLL) resummed DY and Higgs cross sections, matched to NNLO. We find that the inclusion of NLP logarithms is numerically more relevant than increasing the precision to N3LL at LP for these processes. We also perform an analytical and numerical comparison of LP NNLL + NLP LL resummation in soft-collinear effective theory and direct QCD, where we achieve excellent analytical and numerical agreement once the NLP LL terms are included in both formalisms. Our results underline the phenomenological importance of understanding the NLP structure of QCD cross sections.


Introduction
The convergence of perturbative calculations of hadronic cross sections in Quantum Chromodynamics (QCD) is key to understanding the phenomenology of particle collisions at hadron colliders. While the perturbative approach relying on the smallness of QCD coupling constant α s is often successful, its value in certain kinematic limits is conditional on the behaviour of logarithmic corrections. A well-known example is the presence of logarithms of the threshold variable ξ in these corrections. These logarithms grow large as ξ → 0, i.e. approaching the threshold boundary of phase space. One may generically express the structure of the dimensionless partonic coefficient function near threshold as ∆(ξ) = The first term on the right is the class of leading-power (LP) contributions, while the second consists of constants localized at ξ = 0. The former contributions are well-studied and originate from the phase-space difference between virtual and real-emission diagrams after the KLN cancellation and mass factorisation of the soft and collinear singularities takes place. While singular in the ξ → 0 limit, they are regularized by the plus prescription, which makes them integrable. Long ago it was realized that the predictive power of the perturbative approach may be saved by resumming the logarithmic enhancements to all orders in α s , first to leading logarithmic (LL) order [1,2], and soon after extended to subleading logarithmic accuracies [3,4]. These initial results have been re-derived and generalized using either different techniques in direct QCD (dQCD) [5][6][7] or the Soft-Collinear-Effective Theory (SCET) approach [8].
The third set of terms in eq. (1) is down by a single power in ξ with respect to the first. These are the so-called next-to-leading power (NLP) contributions and are the main focus of this work. On the formal side, conjectures on the exponential structure of NLP contributions have already appeared in various places [9][10][11][12][13][14][15][16][17], and have been been proven at NLP LL accuracy for the DY and Higgs production cross sections using two methods. In ref. [18] NLP LL resummation was achieved for general colour-singlet production processes using diagrammatic techniques in the dQCD framework. NLP resummation is also achieved in SCET, both for the DY [19,20] and Higgs [21] production processes 1 . Note that these methods resum NLP LL terms that originate from next-to-soft gluon emissions only. A general resummation framework for soft-quark emissions is still under development (see e.g. [30][31][32] for recent process in the SCET framework).
Despite this interest in the behaviour of the NLP terms, relatively few studies focus on their numerical/phenomenological aspects. Early numerical studies on the threshold approximation of the Drell-Yan (DY) production process already noted the importance of subleading-power logarithms [33,34]. Threshold expansions of the Higgs next-to-next-to-leading order (NNLO) and N 3 LO coefficient functions are studied in [35][36][37], showing that the convergence of this series is slow, thereby stressing the importance of NLP contributions. In ref. [9], the LP resummation formalism was extended to include subleading contributions for the DY and Higgs production processes, and found that these impact the final result significantly. In ref. [38][39][40][41] it was found that the way one handles NLP contributions results in ambiguity of the LP resummed result, and that the numerical difference that follows from this ambiguity is of phenomenological relevance for both DY and Higgs 2 . In [42,43] numerical effects of subleading powers were studied for prompt photon production. All these studies point in the same direction: although not as divergent as the LP contributions, numerically NLP terms can become quite important.
Motivated by these considerations, we aim to answer three questions in this paper. First, what is the behaviour of the threshold expansion of the DY and Higgs partonic coefficient functions? Second, what is the numerical impact of NLP LL resummation in dQCD, and how does it compare to that of subleading logarithmic corrections at LP? And finally, to what extent do the dQCD and SCET NLP resummation formalisms agree, analytically and numerically?
The outline of this paper is then as follows. In Section 2 we study the convergence of the fixedorder DY and Higgs expansions at threshold; for the DY threshold expansion this, to the best of our knowledge, has not been shown earlier (neither for the diagonal nor off-diagonal production channels). For the Higgs case in the dominant production channel, we confirm results that have appeared earlier [35,44]. In Section 3 we review and slightly improve LP and NLP resummation for these processes, and derive the resummation exponent (eq. (49)) that is used for our numerical studies. These are presented in Section 4, where we show the impact of the NLP LL resummation on top of the LP NNLL( ) resummed and matched DY and Higgs cross sections. We also present an assessment of the numerical importance for off-diagonal emissions, for which no resummation has been achieved yet 3 . In Section 5 we perform a careful comparison of dQCD and SCET resummation to LP NNLL plus NLP LL accuracy. We summarize and conclude in Section 6, while certain technical aspects and collections of relevant perturbative coefficients can be found in three appendices.

Threshold expansions at fixed order
In this section we study the convergence of the threshold expansion in fixed order calculations for the single Higgs and the DY processes, first at the parton level, and subsequently for the hadronic cross section. The partonic flux plays an important role in the threshold expansion of the latter. While earlier studies (see e.g. [9,[33][34][35][36][37][38][45][46][47][48]) have addressed partonic and/or hadronic threshold expansions for one or both of these processes, we provide here an extended analysis tailored to study NLP corrections. As such it also serves as an introduction to the numerical study of NLP effects in resummed cross sections in section 4.

Threshold behaviour of the partonic coefficients for single Higgs production
We start with the study of the behaviour of the partonic coefficient functions for single Higgs production at NLO [49][50][51] and NNLO [47,52,53]. First we introduce some relevant quantities and variables. The hadronic cross section for this process is given by 2 We comment on the comparison with our work in section 4.4. 3 Although important progress has been made, see e.g. [16,30,32] where τ = m 2 h /S, with m h the Higgs boson mass and S the hadronic center-of-mass (CM) energy squared, f i (x, µ) are the parton distribution functions (PDFs) and ∆ h,ij (z, m 2 h /µ 2 ) is the partonic coefficient function. Note that the leading order (LO) partonic cross section σ h 0 is factored out. In the infinite top-mass limit 4 it reads (see e.g. ref. [52]) in units of GeV −2 . The strong coupling is denoted by α s ≡ α s (µ), N c is the number of colours, and Fermi's constant G F has the value 1.16639 · 10 −5 GeV −2 . Throughout this paper we employ the central member of the PDF4LHC15 NNLO 100 PDF set [54] for proton-proton collisions with √ S = 13 TeV, corresponding to α s (M 2 Z ) = 0.118 with the Z-boson mass M Z = 91.18 GeV. In practice we use a polynomial fit of these PDFs (see appendix A), which is particularly convenient for obtaining the resummed results presented in section 4. We choose the renormalisation and factorisation scale equal and denote both as µ. The threshold variable is 1 − z = 1 − m 2 h /s, with √ s the partonic CM energy. The hadronic cross section can also be expressed as in terms of the parton flux and where we have suppressed scale dependence in both the flux and coefficient functions. The perturbative expansion of the partonic coefficient functions ∆ h,ij (z) reads h,ij (z) + ∆ (1) h,ij (z) + ∆ (2) h,ij (z) + . . . , where ∆ (n) h,ij (z) includes α n s . The definition in eq. (4) implies that ∆ h,ij (z) is normalized as since the Born-level process involves the gluon-gluon fusion channel. Near z = 1 the function ∆ h,ij (z) can be expanded as The c LP nm are the coefficients of the LP contributions (as in eq. (1)), which contain a factor of (1−z) −1 , while NLP contributions at O((1−z) 0 ) have coefficients c NLP nm ; in general N i LP contributions contain a factor (1 − z) i−1 . Explicit forms of the NLO and NNLO coefficients ∆ (1) h,ij and ∆ (2) h,ij are obtained from the functions η (n) ij (z) eq. (45)- (57) in ref. [52]. We note that the factor of 1/z appearing in eq. (15) of that reference is included in our partonic coefficient functions, i.e. we have 5 In our threshold expansion the factor 1/z is expanded too, an approach also taken in [47], following [48]. One may choose to keep that factor unexpanded. This would lead to somewhat different results, a consequence of truncating the expansion eq. (8). Let us comment here further on the role of the extra factor of 1/z. In general, one may define the partonic coefficient functions to contain additional powers of 1/z, provided a compensating factor is introduced. For some power p we would find instead If one expands the term in square brackets around threshold and truncate that series at some power in (1 − z), the above result is no longer equal to eq. (4), and is thus sensitive to the p additional inverse powers of z that are expanded This makes results of threshold expansions somewhat ambiguous, an aspect also discussed at length in ref. [44] (see in particular the caption of fig. 2 in that reference). In our case we fix the power of 1/z by requiring the universality of LL NLP terms in the dominant channel of colour singlet production processes, as shown in [55,56] at fixed order, as well as in [18] in a resummation context. In both DY and Higgs production the coefficients of the highest power of ln(1 − z) at LP and NLP are, at each order in α s , identical but with opposite sign. This follows from both terms resulting from multiplying the same residual collinear singularity with an overall -dependent power of (1−z). This singularity, which is absorbed in the PDF by mass factorisation, is associated to the lowestorder splitting kernel, such that the logarithms in the finite part are indirectly generated by the splitting kernel too, as is well known. Expansion of the lowest order splitting function, P qq and P (1) gg for DY and Higgs respectively, reveals that they indeed obey this relation between the LP and NLP terms The definition of the (partonic) cross section in eq. (4) satisfies the condition, as we show explicitly for the NLO coefficient function ∆ (a) NLO gg-channel. Diff. Diff.
(d) NNLO qg-channel. The bottom panes show the difference between the truncated and unexpanded expressions. The yaxis range of these difference plots is fixed at 20% of the range of ∆ (1) h shown in the upper panes, for a uniformized comparison. We use eq. (45), (46) and eq. (48)-(50) in ref. [52]) for the unexpanded form of the partonic coefficients. Note that we have included the (trivial) LP approximation of the qg-channels, which vanishes, as these channels start contributing at NLP.
where we have expanded around threshold in the second line, and where we observe the advertised relation between the coefficient of the LL at LP and NLP. We now examine how well the partonic coefficient functions are approximated by the threshold expansion at NLO and NNLO order, for each partonic channel. Our default scale choice is µ = m h . In fig. 1a and b we show the threshold expansion of the two NLO coefficients ∆ (1) h,gg and ∆ (1) h,qg . The third partonic NLO channel (qq) only contributes at N 4 LP, hence we leave it out of our discussion 6 .
We show the results without the δ(1 − z) term, which is present in the gg-channel. In fig. 1c and d we show the equivalent NNLO results. Considering the gg-initiated channel first ( fig. 1a  and fig. 1c), we observe that the LP threshold expansion of ∆ (1) h,gg deviates considerably from the unexpanded result in the z → 1 limit, which might be surprising at first sight. This is caused by the NLP terms: ln i (1 − z) terms do not vanish as z → 1, and are not captured by the LP truncation of the matrix element. Although subdominant to the LP contribution, they are not altogether negligible as z → 1. The unexpanded NLO result is however well-described by the NLP approximation for z 0.2. None of the truncations captures the behaviour below z 0.2 well, due to the factor of 1/z in the partonic coefficient function. The threshold expansion of the NNLO coefficient function ∆ (2) h,gg performs worse ( fig. 1c) than that of the NLO coefficient function: the LP approximation underestimates the unexpanded result in the large-z domain, while the NLP approximation overestimates it. Convergence is seemingly only obtained at NNLP.
Turning to the qg-channel ( fig. 1b), whose LP approximation vanishes, we see that the NLP approximation overestimates the unexpanded result in the large-z domain. Also here the small-z domain is poorly described by the z → 1 expansion of the full partonic coefficient function due to a 1/z factor. The NNLO qg coefficient function ( fig. 1d) shows similar behaviour.

Threshold behaviour of the partonic coefficient functions for Drell-Yan
Next we perform the same studies for the Drell-Yan process. The distribution in the squared invariant mass Q 2 is given by with z = Q 2 s . The α s expansion of ∆ DY,ij (z) is similar to eq. (6), and the threshold expansion of ∆ DY,ij (z) is as in eq. (8). The coefficient ∆ (0) DY,ij (z) is given by δ iq δ jq δ(1 − z) by extracting in eq. (14) the prefactor (see appendix B) where α EM = 1/127.94 is the electromagnetic fine-structure constant at the scale M Z . To compare directly with single Higgs production we set µ = Q = m h for the discussion in this subsection. In fig. 2a and fig. 2c, we show the threshold expansions truncated up to NNLP of the NLO and NNLO qq-channel, whose exact expressions are obtained from ref. [57]. Both perturbative orders are well described by the NLP approximation for a large range of z-values. The convergence of the threshold expansions of the qg-channel, shown in fig. 2b and fig. 2d, is worse than for the dominant qq-channel. In contrast to the Higgs case, the NLP approximation underestimates the ∆ (1) DY,qg and ∆ (2) DY,qg coefficients near z = 1. The NNLP expansion behaves slightly better, although convergence is again slow for this channel. As for the Higgs case, the partonic channels (qq, gg) that open up at NNLO do not contribute at NLP, as these channels require the emission of two soft quarks.
Recapitulating, for both single Higgs production and DY we have seen that the partonic threshold expansion works best for the dominant production channel. The NLP truncation overestimates to a factor (1 − z) 2 from the matrix element. An extra factor of (1 − z) follows from the phase space measure. Diff.
(b) NLO qg-channel. (c) NNLO qq-channel.  the coefficient functions for intermediate and large values of z in these dominant channels, but convergence is reached at NNLP. In the qg-channels, the NLP expansion overestimates the exact result for Higgs production for z → 1, while it underestimates it for DY. At NNLP no substantial improvement is obtained in this channel. The contribution in the z → 0 region, corresponding to the high-energy limit, is for Higgs production more pronounced than for DY, for both production channels. This is due to the factor of 1/z that is part of the partonic coefficient function of the former process, which magnifies small-z contributions. To what extent the various differences also manifest themselves in hadronic cross sections depends of course on the parton flux. This question is addressed in the next section.

Threshold expansions of the Higgs and DY hadronic cross sections
In the previous subsection it became clear that the Higgs partonic coefficient functions have a more pronounced small-z contribution than those of DY. In the hadronic cross sections, these coefficient functions are weighted by the parton flux as in eqs. (4) and (14), which possibly affects the quality of the threshold expansion for the hadronic cross section.
For non-singular terms in the partonic coefficient functions (those that do not contain plusdistributions), the convolutions in eq. (4) and eq. (14) correspond to a point-by-point multiplication  with the parton luminosity function. The product, which is the integrand for the inclusive cross section at fixed τ , is a function of z only. This weight distribution shows which parts of the coefficient functions are enhanced or suppressed by the parton flux, and thus gives much information about how well the threshold expansion can approximate the exact result. We emphasize that for understanding the quality of the threshold expansion, only the shape of the weight distributions for the non-singular contributions matters; the plus-distributions cannot be expanded any further in the threshold limit. For completeness, we nevertheless review the weight distributions for these LP terms briefly, as these notably involve the derivative of the parton flux, as well as the separately integrated partonic coefficient function (see appendix C). For brevity we denote the latter by In fig. 3 we plot the resulting weight distribution for DY and Higgs production, up to NNLO and for Q = µ = m h . To aid comparison we normalize the plotted functions (generically denoted by f (z)) as such that the absolute area under each curve equals unity. The qq parton luminosity is defined as a charge weighted sum of symmetrized parton flux contributions from the five lightest quark flavours: with e q the fractional charge of the quark, normalized to the electromagnetic charge e. The differentiated qq-flux highlights the small-z region more than the differentiated gg-flux, as shown by the dashed lines, but since the LP terms of the partonic coefficient functions for DY are small in that regime (see fig. 2a and fig. 2c), this affects their product with the coefficient functions (solid lines) only little. The product diverges (integrably so) for both channels for z → 1, as one would expect   for the LP terms. Overall, the LP behaviour of both processes is very similar. Around z ∼ 0.8 a small difference appears, caused by the stronger decrease of the derivative of the qq-flux. This results in a slightly more spread-out (i.e. less threshold-centered) weight distribution for DY, but as we mentioned above, this cannot affect the quality of the threshold expansion. We point out that, somewhat counter-intuitively, the LP terms make up only a modest fraction of the total cross section, for both processes (we shall quantify this in fig. 8 and fig. 9).
In fig. 4 we show the normalized non-singular part of the partonic coefficient functions up to NNLO, weighted by the respective parton luminosity (solid lines) 7 . We focus on the dominant production channel for each process and show the z-dependence of the parton fluxes themselves as well (dashed lines). In addition to the natural comparison scale for on-shell Higgs production of Q = µ = 125 GeV ( fig. 4a), we also consider one off-shell Higgs-production scenario with Q = µ = 500 GeV ( fig. 4b). In both cases we see that the product of the parton luminosities and the coefficient functions are suppressed in the small z regime, as a result of the small parton flux. This, in turn, is due to the momentum fractions x 1 and x 2 being large when z → τ since x 1 x 2 = τ z . The individual sea quark and gluon PDFs are small in the large x domain, suppressing the luminosity 8 . Larger values of µ in the flux strengthen this suppression. Therefore, the singular behaviour of the coefficient functions for single Higgs production near z = 0 is more suppressed by the parton luminosity as µ increases. Thus at Q = µ = 125 GeV a notable peak is present in the small z region (it is more pronounced at smaller µ values). It disappears due to the above-mentioned suppression for larger µ values, as seen in fig. 4b for Q = 500 GeV. However, at Q = 125 GeV this feature should still affect the quality of the threshold expansion for the hadronic Higgs cross section.  As the DY partonic coefficient function does not show singular behaviour in the small z-regime, the µ-dependent suppression of the parton luminosity has less impact and we expect the quality of the threshold expansion to be mostly Q-independent. Based on these considerations we therefore expect that, especially for small-Q values, the quality of the threshold expansion of the dominant channel for DY is better than that for Higgs production.
To test this expectation, we compare the parton-luminosity-weighted NLP and NNLP approximations of the partonic coefficient functions to the unexpanded NLO+NNLO result. This product, being the (approximated) integrand of the non-singular (NS) contribution to the hadronic cross section, is denoted by Fig. 5 shows normalized plots of these quantities. We see that the weight distribution for DY production is approximated well, for both Q values, by the NLP truncation over the full range of z.
At NNLP the agreement with the exact result is excellent. For Higgs production at Q = 125 GeV, we again observe that the expansion does not capture the small-z region well, neither at NLP nor at NNLP. The stronger parton luminosity suppression at Q = 500 GeV does aid the convergence, with a NNLP truncation that is almost as good as for DY. Therefore, these plots confirm the expectation that the threshold expansions works better for the dominant production channel of DY than for Higgs production. In fig. 6 and fig. 7 we show similar plots for the qg-channels of both processes. In particular, fig. 6 shows that the z-dependence of the qg-flux is not qualitatively different from the qq and gg flux 9 . Just as for the gg-channel, the parton luminosity suppresses the small-z domain, such that   the enhancement in that region from the factor of 1/z in the Higgs partonic coefficient function is tempered. From fig. 7 we see that the power expansions approximate the full result less well for these channels than for the dominant production channels. Moreover, while the NLP and NNLP approximations for the Higgs qg-channel get significantly better for higher Q values, those for the DY qg-channel do not, and the threshold expansion for Higgs outperforms the one for DY at Q = 500 GeV. Again this is due to the parton luminosity suppression in the small z-region. For Higgs production the truncated expression deviates most from the exact (N)NLO coefficient function in this region, as seen in fig. 1, while for DY ( fig. 2) this deviation is more spread out and therefore benefits less from going to large Q values. We note that the results presented in this subsection contain a level of detail that is of course lost upon integration over z. As such, the integration over z may lead to a seemingly contradictory result: cancellations between under-and overestimations of the exact result across the z domain may cause crude approximations to look more favourable than expected. This is what we observe in the next subsection, where we consider the quality of the threshold expansion of the integrated hadronic cross section. For example, for Q = 125 GeV the integrated NLP truncation approximates the exact NLO+NNLO result for the qg-channel in Higgs production better than the NNLP one, as will be seen in fig. 8, since the more severe overshoot for z 0.3 provides a better cancellation for the undershoot at small z.

Convergence of the threshold expansion in integrated hadronic cross sections
We now turn our attention to the behaviour of the threshold expansions of the total hadronic cross section, starting with single Higgs production. We show the expansions up to N p+1 LP, which follows from  with c to c NLP nm , etc. As a possible δ-contribution cannot be expanded around z = 1, we show this term separately if its coefficient is non-zero. As in the previous subsection, we add the NLO and NNLO contributions in eq. (20), and we restrict our discussion to those partonic channels that contribute at NLP.
The truncated cross sections can be seen in fig. 8, where the power expansions are shown up to N 20 LP. On the left-hand side of this figure, we show results for the gg-channel (up to NNLO) for Q = 125 GeV. The LP expansion severely underestimates the unexpanded part of the hadronic cross section (without the δ-contribution), while the NLP expansion overestimates it by a much smaller amount. After including the NNLP term, the expansion stabilizes 10 . We also observe a significant δ-contribution, as is well known for this process.
For the qg-channel (right-hand side of fig. 8), we observe that the NLP truncation exactly produces the unexpanded cross section. This is however somewhat of a coincidence: the (N)NLO NLP truncation underestimates (overestimates) the magnitude of the negative unexpanded (N)NLO contribution by a similar amount. When added together, these under-and overestimates cancel. The NNLP truncation for the integrated cross section performs much worse, consistent with what we predicted based on fig. 7a. An overestimate of the negative contribution from this channel perseveres for higher-power truncations, and the expansion converges only very slowly to the full result. Again, this is due to the NLO and NNLO coefficients having a negative contribution at z → 1, compensated by a large positive contribution for z → 0, which is however only slowly reconstructed in a 1 − z expansion.
Also for DY production at Q = 125 GeV, the LP expansion of dominant channel is not a good approximation of the exact cross section, as seen in the left-hand side of fig. 9. The NLP expansion shows moderate overestimation, as seen in fig. 5a as well, but performs significantly better already. In general the power expansion converges quickly. For the qg-channel (right-hand side of fig. 9) we see that the NLP (NNLP) expansion overestimates (underestimates) the absolute value of the  unexpanded NLO+NNLO coefficient, consistent with the deviations observed in fig. 7a. Contrary to the Higgs qg-channel, the threshold expansion does converge after including the N 4 LP contribution. Before we conclude this section, we comment on the behaviour of the threshold expansions for other values of Q. For the gg-channel in single Higgs production, the convergence of the threshold expansion happens faster for higher values of Q, which is a direct consequence of the stronger 1/z suppression of the gg luminosity function. Also for the Higgs qg-channel the threshold expansion converges more slowly for small Q values than for larger ones, for the same reason. On the other hand, the behaviour of the threshold expansion of the qq and qg DY channels marginally improves at higher Q values. We conclude that the threshold expansion for NNLO cross sections is in general more reliable for DY than for single Higgs production for both the dominant (qq/gg) and subleading qg production channels.

NLP resummation in dQCD
Having exhibited the effect and quality of NLP approximations in fixed order cross sections, we now turn to consider NLP effects for resummed cross sections. This section discusses analytical aspects of NLP resummation in direct QCD (dQCD), whereas numerical results are shown in Section 4. We concentrate on the dominant channels, and address LP and NLP resummation at the same time.

From LP to NLP resummation
Resummation in dQCD is customarily performed in Mellin-moment (N ) space, where the cross section is a product of N -space functions. Thus, for the gluon-fusion contribution to Higgs production in eq. (2) one obtains and similarly for the DY production (with g → q and σ h 0 → σ DY 0 ). We use a boldface notation for Mellin-transformed quantities. To perform the resummation in N -space, one uses the resummed perturbative coefficient ∆ aa (N ). The resummed hadronic cross section in momentum space is obtained after taking the inverse Mellin transform where we choose the Minimal Prescription [58] for the N integration contour. This corresponds to choosing c = C MP to left of the large-N branch-cut (starting for us at the branch-pointN = exp [1/(2α s b 0 )]) originating from the Landau pole, but to the right of all other singularities. We note that this approach works for both LP and NLP terms, as the branch cut is identical for both terms.
To minimize numerical instabilities, the integration contour is usually bent towards the negative real axis. We include also N -space PDFs in our inverse Mellin transform, for which we use a fitted form of the PDFs (see the discussion around eq. (131)). In ref. [18] a resummed partonic coefficient for colour-singlet production processes was derived that also incorporates the LL NLP contributions for the dominant channel The exponent E LP+NLP a contains the diagonal DGLAP splitting function P aa , expanded up to NLP in the threshold variable. The factor of 2 multiplying this term reflects the initial state consisting of identical partons (either gg or qq). The soft wide-angle contributions are collected in D aa (α s ). Both terms are process-independent to the extent that they only depend on the colour structure of the underlying hard-scattering process. The overall +-subscript denotes that the plus-prescription needs to be applied to all LP contributions. After the Mellin transform has been carried out, the D LP+NLP aa and E LP+NLP a terms contain all logarithmic N dependence at LP and leading-logarithmic N dependence at NLP. It also contains some contributions that are beyond LL NLP (i.e. the argument of α s in the function D aa creates NLP contributions beyond LL, as we will see, and so does the 1/z-dependence of the upper limit of the k T integral). The process-dependent function g 0 (α s ) collects the N -independent contributions. For an N k LL resummation, we need this function only up to O(α k−1 s ), but if available we include the O(α k s ) terms as well, which upgrades N k LL to N k LL resummation. The splitting function P LP+NLP aa can be written in the LP+NLP approximation as where 11 The coefficients A (n) a are known up to fourth order [59], but for NNLL accurary we need them only up to third order [60,61]. They are given explicitly in appendix D. The perturbative expansion for the function D aa (α s ) reads which starts at the two-loop level since D aa (see appendix D) is needed for a resummation at NNLL( ) accuracy.
Upon substitution of eq. (24) into eq. (23), we may compare it with its LP counterpart [3,4] ∆ dQCD, LP where the plus-prescription is already applied. We highlight two changes with respect to (23). The first is that the splitting function is approximated to LP accuracy instead, removing the additional −1 term of eq. (24). The second is that the upper limit of the k T integration on the second line, reflecting the exact phase-space constraint on the soft emission, has been replaced by its LP approximation. The same replacement has been made in the argument of α s in D aa on the first line. How to calculate the integrals in eq. (27) is outlined in ref. [4,6,62], and extended to accommodate NLP contributions to the splitting function in ref. [11]. However, the latter reference did not implement the exact phase-space constraint. If one does, a simpler formula can be derived that obtains the NLP contribution directly from the LP exponent using a derivative with respect to the Mellin moment N , which we show in the next subsection.

The NLP exponent through differentiation
Let us calculate E LP+NLP a at LL accuracy, since we are interested in the NLP LL contribution for which we only need A The integral over k T may be evaluated using the QCD β-function, for which we only need the one-loop coefficient b 0 (see appendix D), and write where we use α s ≡ α s (µ 2 ). This may be written as where we have separated the logarithms with scale-dependence from those with z-dependence. Since we wish to resum the highest power of the threshold logarithms at each order in α s , we select the pure ln k ((1 − z) 2 /z) term from the above expression, discarding the explicit scale logarithms that only contribute at NLL accuracy and beyond (we would have obtained the same result by making the replacement µ 2 → Q 2 in the lower integration boundary of eq. (28)). Expanding and retaining the first two terms, Eq. (28) becomes In the second line we introduced shorthand notation for the Mellin integrals, which are evaluated through their generating functions G F . We define for F = {D, J , J } and This method was proposed in [6] for the LP contributions D k and generalized to the NLP integrals J k in [11]. The terms labeled by J k were not included before. Although of sub-leading logarithmic accuracy, they are important for our final result. Expanding eq. (34) around the limit N → ∞ yields, We see that pure LP contributions in z-space, as contained in D k , give NLP contributions in Mellin space (i.e. terms proportional to 1/N ). In dQCD, the resummation is done in Mellin space, and in what follows, when we refer to 'NLP' in dQCD, we mean terms that are proportional to 1/N . We now observe that all O(1/N ) terms can be generated from the LP contributions in N -space by means of a derivative with respect to N Note that this only holds when including J , as it cancels the η 2 term in the square brackets of eq. (35a). Taylor expansion of the LP contribution (i.e. the term in eq. (36) without the derivative) around η = 0 yields with Γ (n) denoting the n-th derivative of the gamma function, i.e.
In this way, we have isolated the η behaviour in one simple factor, such that the derivative of eq. (33) amounts to the relation Rewriting eq. (36) using eq. (37) we obtain, via eq. (33), for eq. (32) Terms in the sum with n > 0 are sub-leading logarithmic terms in N , but some of those are easily included by redefining N toN = exp [γ E ] N . To illustrate this, we note that where all corrections contain constants of higher transcendental weight (ζ(n) for n ≥ 3). Using eq. (41) and shifting the summation index n → n + 2 for the second term in eq. (41) (where n = 0 and n = 1 do not contribute) yields having recognized the binomial series for (ln N + γ E ) k±1 ≡ ln k±1N in the last line. We see that the ζ(2) term in eq. (41) contributes only at NNLL accuracy. Defining λ = b 0 α s lnN and performing the summation over k, this factor results in The +1 contribution in the last line is included in the O(α s ) contribution of g 0 . The first term in the square brackets is instead included in the NNLL contribution to the resummed exponent (see eq. (152) of appendix D). At LL accuracy we have where the LL NLP resummation function h a is obtained from the LL LP function g a by taking the derivative towards N . Let us stress that the logarithmic accuracy of the generated NLP terms is limited by that of the LP function on which it acts. Therefore, this result is strictly of LL accuracy at both LP and NLP. We will use the resulting NLP function for the numerical studies in section 4.
By using the full A(α s ) function for E LP+NLP a the result is straightforwardly extended to higher logarithmic accuracy at LP, and including partial N i LL NLP terms through the derivative in N of the relevant LP resummation functions A somewhat different form of this LP + NLP result proves to be useful in section 5, which is obtained by acting with the derivative on the result of ref. [62] where the Mellin transform has been expressed differently This is valid to arbitrary LP logarithmic accuracy, and to NLP LL accuracy, as we do not have control over subleading logarithmic NLP contributions that originate from wide-angle soft gluon emissions.
For the wide-angle contribution in eq. (23), we can use the same method to arrive at a derivative formula, but where the derivative term changes sign For the first non-zero term in the perturbative expansion of D aa the result is This is an LP NNLL contribution that does not give rise to LL contributions at NLP, hence we may drop the derivative term. Our final form for the dQCD resummation exponent accurate to LP NNLL and NLP LL order reads The wide-angle contribution of eq. (48) is contained in g a . Explicit expressions for these resummation exponents are collected in Appendix D. We will use eq. (49) for our numerical studies in section 4.

Numerical studies of LP and NLP resummation cross sections in dQCD
Having set-up the formalism for dQCD resummation at NNLL at LP and NLP LL, we turn to the numerical study of this formalism for various processes, in the context of LHC collisions at √ S = 13 TeV. We use fitted PDFs that allow for an analytical evaluation of the Mellin transform, as explained in appendix A. For DY and Higgs production we will show resummed observables that are matched to the fixed-order result at NNLO. The matching is defined by Note that we include the full fixed-order result in σ (fixed order) , i.e. also the sub-dominant production channels. The second term on the right is the expanded resummed observable. We calculate this term by Taylor-expanding the resummed coefficient function in N -space where we only need the terms up to O(α n s ) for matching with an N n LO fixed-order calculation. The result is substituted into eq. (21), after which the Mellin-space inversion is handled via eq. (22). Terms that are of higher power than O(1/N ) are created in the expansion (51). They are kept in the matching and thus subtracted from the resummed result as they are also contained in the complete fixed-order expression. We perform the LP resummation at NNLL accuracy, where NNLL resummation offers an improvement over NNLL by the inclusion of the exact N -independent terms at NNLO (i.e. we include g 0 up to O(α 2 s )). This is not strictly necessary for NNLL resummation, but it is relevant in case of large virtual corrections at the two-loop level, such as for single Higgs production. We stress again that we only resum the channels that contribute at LP. To examine the impact of NLP resummation on colour-singlet processes other than Drell-Yan or single Higgs production, we also show at the end of this section results for di-boson and di-Higgs production, which we include to LP NLL accuracy and do not match to fixed higher-order results.

Single Higgs production
To obtain the total resummed Higgs cross section in N -space, we start from eq. (21) with σ h 0 given in eq. (3) and the resummed contribution ∆ dQCD gg (N, Q 2 /µ 2 ) of eq. (49) with a = g. As mentioned above we include g 0 up to O(α 2 s ). We vary Q(= m h ) with the aim of exploring the resummation effects more widely than only for the physical scale Q = 125 GeV.
The results are shown in fig. 10, where all lines have been normalized either to LO (a) or NNLO (b). All resummed results are matched to the NNLO fixed-order result. The large enhancement due to the NNLL matched resummation with respect to the LO result is caused by the sizable δ(1 − z) contribution to both the NLO and NNLO fixed-order cross sections, which enter the resummed distribution via g 0 (α s ). The enhancement of the LP NNLL + NLP LL resummed result, with respect to the LO distribution (using µ = Q), is roughly 300 − 360% in the considered Q range, the strongest increase occurring for smaller Q values. When compared to the matched LP NNLL resummed cross section, we find an NLP enhancement between 4.3 − 6.3%, with larger effects for smaller Q values. This increase can be contrasted with the effect of increasing the LP accuracy to N3LL. To obtain this we need to include only the g a (λ, Q 2 /µ 2 ) function in the resummation exponent 12 , given that g 0 (α s ) is already included up to O(α 2 s ) for NNLL accuracy. Only a modest (negative) N 3 LL correction to the NNLL resummed cross section is found: between −(0.3 − 0.5)%. We have verified that the NLP corrections are competitive with the numerical increase from NLL to NNLL, and we show explicitly that they dominate over the increase from NNLL to N 3 LL 13 . The scale uncertainty increases somewhat after including the NLP resummation. For the NNLL resummed result, the scale uncertainty is between −0.03% and +6.2%, whereas the NLP resummed result shows a scale uncertainty between −2.9% and +8.3% (we expect that the inclusion NLP NLL terms would reduce this scale uncertainty). Below, we report explicitly on the Higgs cross sections at Q = 125 GeV, corresponding with the physical Higgs mass σ (NNLO) pp→h+X (m h ) = 39.80 +9.7% −9.0% pb, 12 We extract this function from the publicly available TROLL code [40]. 13 One could also compare to the N 3 LL results, where one includes the O(α 3 s ) contribution to g0(αs). Here we choose not to do so, in order to directly compare the numerical contribution of terms that are of O(α k s ln k−2 (N )) with those that are of O(α k s ln(N ) k /N ) (where k runs from 1 to ∞). If we would compare at LP N 3 LL instead, the O(α k s ln k−2 (N )) terms in the exponent would be multiplied by a different hard function than the LP NNLL + NLP LL result, which pollutes the comparison with a possibly sizeable constant contribution. We thus find a notable NLP LL contribution in the diagonal channel. Of course, the resummation of the qg-channel is missing, which would also result in NLP LL enhanced terms (see fig. 8), and may potentially alter the observed NLP effect. To estimate the size of such contributions we include the NLP LL contribution from the qg-channel at N 3 LO, the order at which we expect the largest contribution from a potential resummation in that channel. The full N 3 LO results for the Higgs total cross section are available in the infinite-top-mass limit [63], and can be inferred from the iHixs2 code [64]. We find that the qg induced NLP LL contribution reads, in Mellin space, which coincides with the prediction of ref. [16]. This qg-contribution results in a correction of the NLP LL resummation effect of −0.5% (−3%) for small (large) Q values (see fig. 11a). The O(α 3 s ) contribution of the qg-channel thus gives a negligible contribution to the NLP LLs. Adding the terms of eq. (53) does not lead to a noticeable reduction of the scale uncertainty (not shown in the figure). Note that the smallness of the qg-initiated result is not caused by the partonic flux: the qg flux exceeds the gg flux at µ = 125 GeV. Instead, this contribution is relatively small because, in contrast to the NLP LL contribution from the gg-channel, the qg-contribution is not proportional to the sizable higher-order constant terms of the LP gg-channel contained in the O(α s ) and O(α 2 s ) contributions to the hard function g 0 (α s ). A similar hard function is not included for the qg-channel, as cannot use an exponentiated form for the qg contribution, nor know what the hard function in that case would be. We expect that the qg-channel will play a larger role in the DY process, where the g 1 0 and g 2 0 coefficients are comparatively small.

DY invariant mass distribution
For the computation of the N -space resummed DY invariant mass distribution we use with σ DY 0 given in eq. (15), ∆ dQCD qq (N, Q 2 /µ 2 ) in eq. (49) and the sum runs over quarks treated as massless, q ∈ {u, d, c, s, b}. We show the resulting ratio plots in fig. 12, where as in the Higgs case, all resummed results are matched to the NNLO result.
The LP NNLL + NLP LL resummation enhances the LO distribution (using µ = Q) by 15 − 34%, with increasing effect for larger Q values. The NLP contribution provides an increase with respect to the LP NNLL resummed (and matched) distribution of only 0.35 − 1.15%. Larger NLP enhancements are found for very small values of Q, where one moves further away from threshold. As for the Higgs production case, the NLP increase dominates over the N 3 LL effect, which is −(0.1 − 0.3)% as shown in fig. 12b, where the larger deviation is only found for small values of Q (as best seen in fig. 11b). The larger size of the NLP LL with respect to the N3LL is not as pronounced as in the Higgs production case.   The scale uncertainty of the NNLL resummed result lies between −4.4% and +5.4% for the range of Q values shown in fig. 12, while that of the NLP resummed result is between −4.8% and +5.8%. Therefore, by including the NLP LL contribution, we find a modest increase in the scale uncertainty of the resummed result, which would expect to decrease if NLL resummation at NLP were available. Note that for large Q values (Q > 1 TeV), the central value of the NLP correction lies outside the LP resummed uncertainty band.
As in the Higgs production case, the NLP LL resummation is not complete, since the qg-channel is missing (see fig. 9). Using the results of ref. [16], we may obtain the N 3 LO NLP LL contribution stemming from the qg-channel, resulting in As already anticipated, in contrast to what was observed for single Higgs production, we see that the addition of this result to the resummed qq-channel significantly alters the NLP effect ( fig. 11b).
The qg-contribution gives a −54% correction to the NLP effect of the dominant channel for small Q values, while for larger Q values it gives a −44% correction. As for Higgs production, we find that the scale uncertainties do not decrease after adding the qg NLP LL contribution at O(α 3 s ). This is perhaps not surprising: for both the Higgs and DY processes we introduce additional scaledependence via α s in either the NLP LL function h (1) a for the leading channel, or via eq. (54)/(55) for the off-diagonal channel. This scale dependence is not balanced at LL, but could (for both channels) be balanced at NLP NLL via the introduction a scale term proportional to ln(Q 2 /µ 2 ). However, since we have no control over any other contribution that might appear at NLP NLL, we refrain from adding this contribution.

Other colour-singlet production processes
To see how general our results for the impact of NLP effects are we extend the analysis to the di-Higgs and di-boson production processes. As it is not our aim to provide a precise phenomenological prediction, but rather analyse the numerical effects of NLP LL resummation, we consider the (unmatched) resummation of these processes up to LP NLL + NLP LL, and include only the leading-order contribution to g 0 (α s ). We start by considering the di-Higgs production process, which is also dominated by gluon fusion. Threshold resummation has been achieved first up to NLL in the SCET framework in the heavy-top-mass limit, including top-quark-mass dependent form factors to partially correct for this [65] approximation. This study has been extended to NNLL in ref. [66], and the inclusion of the full top-mass effects was studied up to NLL+NLO in ref. [67]. Using full top-mass dependence, which we study here, the lowest-order expression for the hadronic differential distribution dσ pp→hh dQ is given by [68] dσ pp→hh with τ = Q 2 /S and where the integration variable t is given by while the integration limits are obtained by setting cos θ = ±1. The exact expressions of ref. [68] are used for C ∆ , C , F ∆ , F and G , where no approximation on the mass ratio between the top-quark and the Higgs boson has been applied. For the resummation of this process, it suffices to follow the same procedure as before. That is, we Mellin transform eq. (56) with respect to the hadronic threshold variable τ . Then we replace the partonic LO coefficient σ gg→hh (Q 2 ) by its resummed version, σ gg→hh (Q 2 )∆ dQCD gg N, Q 2 /µ 2 . This replacement is valid at LP (see ref. [69]). We exploit the universal structure of next-to-soft gluon emissions [55] to justify the same substitution at NLP. The resummed expression then reads Since we work at NLL accuracy, we use the expression of eq. (49) for ∆ dQCD gg N, Q 2 /µ 2 with g 0 (α s ) = 1 and neglect the g (3) contribution in the exponent.
The results are shown in fig. 13, where we used a factorisation/renormalisation scale µ = Q/2, which is the scale for which higher-order corrections are smallest [70]. We see that the LP NLL result gives a correction to the LO distribution of −3.2% to 4.7%, while the NLP LL terms cancel the partially negative LP correction, leading to a substantial enhancement of the LO distribution of 12.6 − 17.8%. Larger corrections are found for higher values of Q. An increase of the NLP LL + LP NLL result with respect to the LP NLL result is found to be between 12.5% and 16.3%. As for the DY and single Higgs processes, we thus again find that the NLP LL effect is sizeable. For comparison, the increase of going from LP LL to LP NLL (not shown here) is of the same order 14 .
We now move on to W + W − and ZZ production. Threshold resummation up to NNLL was considered in the direct-QCD framework in ref. [71], and for the SCET framework in ref. [72,73]. Transverse-momentum resummation for both processes has also been performed in ref. [74]. We use the LO expressions from ref. [75,76] 15 . We may write, similarly to eq. (56) with the LO partonic cross sections given by where The expressions for the coefficients c tt i , c st i , and c ss i , and the functions F i , J i and K i may be found in ref. [75]. Note that the partonic coefficient functions σ q iqi →V V depend on the left-and right-handed couplings of the quarks with the Z-boson, which are of course different for up-type quarks and down-type quarks. Following the same procedure as for the di-Higgs results, we obtain a resummed expression for both cases that reads The results are shown in fig. 14. The two processes are identical from a resummation point-ofview, therefore, we find an LP NLL increase of the LO distribution between 12.2 − 20.8% for both processes, whereas the NLP LL + LP NLL increases the LO distribution by 19.2 − 26.2%. As for the other processes, larger enhancements are found for higher values of Q. The increase induced by the NLP LL contributions with respect to the LP NLL result is between 4.5 − 6.8%, which is smaller than the correction that was found for the di-Higgs process. This is not surprising: it is well known that gluon-initiated processes receive larger LL threshold corrections due to C A > C F . The difference between LP LL and LP NLL resummation for the V V processes is around 1%, which is not dissimilar from the increase found for the inclusion of the NLP terms. We therefore again see the numerical importance of the NLP LL terms.
What we have observed for all processes is that the NLP LL contribution is not negligible, and in many cases exceeds the size of higher-logarithmic LP contributions. To understand this, note that the saddle-point approximation for the inverse Mellin transform integral leads to values ofN around 1.5 − 3 [38,39,77,78]. With α s 0.12, we have α s ln(N ) α s ln(N )/N , while the NNLL term α 2 s ln(N ) is suppressed by α s . Naturally, we do not suggest to cease pursuing effects of higherlogarithmic LP contributions. Firstly, one can include the O(α s ) contribution in the hard function g 0 (α s ) (i.e. leading to LP NLL ) accuracy), resulting in a large correction if the constant contribution of the NLO correction is sizeable, as we have seen for single Higgs production, and also seems to be the case for di-Higgs production (see e.g. [67,79,80]). Secondly, although the central result might change minimally by including higher-logarithmic contributions at LP, such contributions stabilize the scale dependence of the prediction.

Some final remarks on NLP LL resummation in dQCD
In this section we have explored the numerical behaviour of the NLP h (1) a correction (eq. (49)) on a few selected colour-singlet production processes. For all processes we reach a similar conclusion: the inclusion of NLP LL resummation leads to a numerically noticeable result. We briefly comment on the relation between the work on NLP resummation presented here and in refs. [40,81] for single Higgs production. The methods employed in this paper may best be compared to their ψ − soft 2 prescription without exponentiation of the constant contributions 16 . The ψ − soft 2 prescription consists of replacing λ in the resummation exponents g (i) (λ) by the combination of polygamma functions λ = α s b 0 (2ψ 0 (N ) − 3ψ 0 (N + 1) + 2ψ 0 (N + 2)). Denoting their resummation exponent at LP NLL with ∆ ψ−soft 2 aa , the difference with our LP NLL + NLP LL resummation exponent (eq. (49) without g (3) a and with g 0 (α s ) = 1) may be evaluated at O(α s ), where we obtain The difference at LP is of NNLL accuracy, and originates from using N exponentation in the ψ − soft 2 prescription rather thenN exponentiation, which is what we have used in our work. At NLP, the difference is of NLL accuracy, whereas the discrepancy at NNLP starts at LL with the term proportional to ln(N )/N 2 . They observe that the correction on the fixed-order single Higgs production cross section is increased with respect to standard resummation by making the ψ − soft 2 replacement, consistent with our observations on the effect of NLP LL resummation. No soft-quark correction was considered in their work. To the best of our knowledge, a similar analysis at the hadronic level for DY has not been performed. Returning to our results, we found in this section that for single Higgs-production the central value of the LP NNLL + NLP LL resummed coefficient lies outside the uncertainty band for the LP resummed (and matched) result. For DY production processes this is observed for Q > 1 TeV. Both processes show that the numerical effect of resumming the NLP LL terms exceeds that of improving LP resummation to N 3 LL accuracy. Scale uncertainties seem to slightly increase after the NLP LL result is included. We expect that the scale uncertainty can be reduced only after the inclusion of the (unavailable) NLP NLL contributions. For di-boson production the numerical enhancement of upgrading a LP LL calculation to NLP LL rivals (in the case of di-Higgs) or exceeds (for W + W − /ZZ) that of going from LP LL to LP NLL, for central scale choices. In general, NLP LL effects originating from next-to-soft-gluon emissions are found to be larger for gg-induced processes than for qq-induced processes. On the other hand, we find that the off-diagonal qg NLP LL contribution is more important for the qq-induced DY process than for the gg-induced Higgs process.

Comparison of NLP resummation in SCET and dQCD
NLP corrections have been calculated and resummed not only in direct QCD, as described in section 3.1, but also within soft-collinear effective field theory (SCET). While the two approaches should formally give the same result, numerical differences have turned out to be sizable for LP resummation, and originate from power-suppressed contributions that are included differently in the two methods [39,81,[84][85][86][87]. Now that the exact LL NLP contribution (in z-space) has been determined in both approaches, it is important to perform such a comparison again, and we do so in this section. We shall see that the inclusion of the LL NLP contribution reduces the numerical differences between the two approaches. For the comparison we consider Drell-Yan and Higgs production, at LP NNLL (i.e. we use eq. (27) with g 0 up to O(α s ) for dQCD) and NLP LL accuracy. We refrain from matching to fixed-order calculations as this has no bearing on the comparison of resummed calculations. Below we first establish notation for SCET resummation, and briefly discuss the results of [19][20][21], revisiting in particular the so-called kinematic correction. We obtain a slightly modified SCET resummation formula, which will ease our task of relating the SCET formalism to dQCD. We then proceed to the comparison of the NLP LL contributions within the two theories. Note that while presenting the analytical comparisons, we use the β-function coefficients b n and β n = (4π) n+1 b n (eq. (146) and (149)) interchangeably, as the b n are often used in dQCD, whereas β n are usually preferred in SCET. The section is concluded by discussing numerical results for both methods.

NLP Resummation in SCET
Using the results of ref. [19][20][21], we write the analogue of eq. (27) in SCET up to NLP in z-space as where ∆ SCET,LP/NLP (z) resums large logarithms at LP and NLP, respectively, and the scale dependence of these factors is suppressed. The NLP term ∆ SCET,NLP (z) consists of two contributions: a kinematic factor, which originates from taking the partonic cross section to LP in the SCET operator product expansion, but retaining the exact kinematics up to NLP, and a dynamical correction that arises by considering sub-leading-power operators We first focus on the kinematic correction, which together with its LP counterpart, is expressed in the factorised form [19,20] ∆ DY,SCET,LP+NLP for the DY partonic coefficient function of the qq-channel, and for the gg-channel contribution to Higgs production, [21]. The hard function H depends on shortdistance physics, and is different for DY and Higgs. The function S 0 denotes the Fourier transform of the position-space soft function where the latter is defined as a vacuum-expectation value of time-ordered products of Wilson lines associated to the directions of the external partons [19,88] S Here C = N c for Drell Yan, and C = N 2 c − 1 for Higgs production. Furthermore, the soft Wilson lines are defined in terms of soft gluon fields in SCET where the light-like vectors n ± are oriented along the directions of the incoming partons. At NLO the Drell-Yan and Higgs soft functions are related by the simple replacement C F → C A . The normalizations for the Higgs and DY case are chosen, as explained in section 2 and further discussed in appendix B, such that we have a universal NLP LL correction for both processes, resulting in a 1/z difference between the Higgs and DY factorisation formulae, as it is evident by comparing eq. (68) with eq. (69). The variable Ω * is the reciprocal energy, and is given by Eqs. (68) and (69) are expanded in powers of 1 − z, which implies in particular that Ω * is expanded around Ω ≡ Q(1 − z). Then, the LP contribution reads where we have suppressed the DY/h labels of H and S 0 . At fixed-order, both H and S 0 depend on the factorisation scale. The hard function contains factors of ln(Q 2 /µ 2 ), and the soft function contains ln(Q 2 (1 − z) 2 /µ 2 ). Therefore, depending on the choice of µ, either H or S 0 will result in a large logarithmic correction. However, the hard and soft functions obey renormalisation-group equations (RGEs). According to the effective-field theory approach, it is then possible to choose a natural scale, µ h Q for the hard function, and µ s Q(1 − z) for the soft function, such that fixed-order logarithms are small. The original large logarithms are then resummed using the RGEs to evolve the hard and soft function from their natural scale to the scale at which the cross section is evaluated. We will recall the resummed partonic cross section one obtains starting from eq. (73) in what follows, but let us first revisit the NLP terms originating from the expansion eqs. (68) and (69), i.e. the kinematic correction.
To obtain eq. (73), we made the choice to expand eq. (68) and eq. (69) around Ω * = Ω and set s = Q 2 . At NLP, this choice leads to the kinematic NLP corrections K 1 , . . . , K 4 where The K 1 term comes from identifying ∂ ∂Ω * = ∂ ∂Ω + O(1 − z), and K 2 originates from the Taylor expansion of the soft function around Ω * = Ω. The K 3 term is the first order of the expansion of Q/z in the case of DY, and Q/z 2 in the case Higgs, hence the resulting contribution S K 3 is different for these processes. Finally, K 4 contains the Taylor expansion of the hard function around s = Q 2 /z. Evaluating these terms up to NLP yields for their sum 17 [19][20][21] S DY K (Ω) = 2 With these definitions, the LP partonic cross section with full kinematic dependence is expanded as where ∆ LP (z) is strictly LP, as defined in eq. (73), and the corresponding NLP kinematic correction is given eq. (74). While this choice is consistent with a systematic power expansion, it is also possible to consider other expansion schemes, for which parts of the kinematic correction are not expanded, but instead kept within the LP term. Here we discuss one such case, motivated by the comparison with dQCD resummation. In this expansion scheme we keep factors of Q/ √ z exact, i.e. we expand Ω * aroundΩ = Q √ z (1 − z). Comparing with eq. (28), we immediately see that this choice consistently keeps the same kinematic factor in the exponent, both in the dQCD and in the SCET approach. The LP factorisation then reads This choice leads for the K 1 and K 4 contributions of the kinematic correction (eq. (74)) to a replacement of Ω →Ω. The kinematic corrections K 2 and K 3 change and becomê (Ω). At the same time, the NLP LL contribution, which is absent for the DY case, and equal to the scale logarithm in the Higgs case, is unchanged. Explicitly, we havê The NLP NLL contribution that was present in eq. (76) is now removed by the new expansion scheme. This contribution represents a phase space correction, which is universal. Instead of residing in the kinematic correction, it is now contained in the modified LP term of eq. (78); in this way it too can be readily resummed. With this result in mind, we move to the resummed SCET LP expression, obtained by evolving the hard and soft function to an equal scale. If we use the LP factorisation as in eq. (73) one obtains the resummed partonic coefficient function while using the modified scheme of eq. (78) we obtain .
(82) In these equations the function η a is fixed to after the derivative with respect to η a has been taken. We provide the process-dependent hard function H(Q, µ h ) and the Laplace transform of the soft function,s a ln Q 2 µ 2 s + ∂ ∂ηa , µ s , for Drell-Yan and Higgs production in appendix D, up to NLO, which is sufficient to achieve resummation at NNLL accuracy. The corresponding expression up to NNLO, necessary to achieve resummation at N 3 LL can be found in [89] and [84], respectively 18 . In eqs. (81) and (82) the actual resummation is organized in the exponent U , given by [89] with U DY ≡ U qq (Q, µ h , µ s , µ) and γ V /S = γ V (see appendix D) for DY production . For the Higgs production case we have γ V /S = γ S and include an additional contribution [84] whose role will become clear later. The resummation factor U aa (Q, µ h , µ s , µ) is expressed in terms of the Sudakov exponent S a and the anomalous dimension exponents a γ , defined as In appendix D we provide the explicit expansion in α s of both S a (Q 2 , µ 2 s ) and a γ i (µ 2 , ν 2 ), as well as the anomalous dimensions they depend upon. Note that the functions A a (α s ) (from dQCD) and as can be seen by comparing the equations in appendix D. Comparing eq. (81) with eq. (82) we see that the difference between the two schemes amounts to a factor of z −ηa , which resums NLP NLL contribution of kinematic origin. The resummation formula defined in eq. (81) coincides with the one found in [84,89]. What is new from our discussion above, is the fact that the NLP LL resummed contribution can be added consistently to both eq. (81), as was done in [19,21], or to eq. (82), as we are doing here. The latter is particularly useful for the comparison with dQCD. The ultimate reason that allows one to proceed in this way is the fact that the NLP LL contribution to the kinematic soft function is unaltered by the different convention that leads from eq. (76) to eq. (80), which affects only NLLs at NLP. Thus, the LL contribution from the NLP kinematic soft function is added unaltered to the LL contribution from the NLP dynamical soft function, according to eq. (67), preserving the universality of the total NLP LL contribution, whose resummed expression takes the form [19,21] where the function S LL (µ 2 , ν 2 ) reads Interestingly, we can directly obtain the NLP LL SCET contribution from the LP LL one. The LP LL SCET contribution is obtained after setting z −ηa = 1,s a = 1 and H(Q, µ h ) = 1. Furthermore a γ = 0 at NLL, such that the evolution exponent eq. (84) simply becomes The NLP contribution can then be obtained from this evolution exponent via a derivative By now evaluating β α s (µ 2 s ) to its lowest order and recognizing that A (1) a = C a , we readily find that this result is equal to eq. (89). It is remarkable that we could perform a similar derivative trick in dQCD (eq. (45)) to obtain the LL NLP results. There the derivative is taken with respect to N , while here it is taken with respect to the coupling evaluated at the soft scale. We will see in what follows that these forms are actually equal.

Analytical comparison of dQCD and SCET at NLP
We are now in a position to relate resummation in dQCD vs SCET, where we closely follow ref. [39,81,85], where this was performed at LP, up to NNLL. We keep the same LP logarithmic accuracy, but we extend their analysis to NLP LL (in N -space). To compare resummation in dQCD vs SCET, it is convenient to convert the SCET resummation formula to Mellin space. We set the hard and factorisation scale (µ h and µ) equal to Q, and to obtain a function whose Mellin transform may be performed analytically, we assume that the soft scale µ s is independent of z [86]. Furthermore, we will work with the expansion scheme in eq. (82), and comment on the role of z −ηa at the end of this section. Then we can write 19 Up to NNLL accuracy we may usẽ Furthermore, the ratio of Γ-functions in eq. (93) may be expanded in the large-N limit, resulting in By using eq. (94) and eq. (95) in eq. (93), we obtain Recognizings a ln 2 Q 2 N 2 µ 2 s , µ s in the second line we can write this as For the NLP SCET result in Mellin space we use eq. (89) and take µ s again independent of z, which results in We now compare eq. (97) and eq. (98) to ∆ dQCD as defined in eq. (23) with E LP+NLP a in eq. (46) and D LP+NLP aa in eq. (47), starting with the former. By defining the derivative operator eq. (46) may be rewritten as where the order of integration of the two integrals is interchanged on the second line. We now split up the integration domain of the two integrals on the second line by introducing an auxiliary scale µ s for the first, and an intermediate scale Q for the second integral. We obtain with where we recognize the SCET Sudakov exponent S a (Q, µ s ) of eq. (86), and We call the latter function E Landau a as it generates the branch cut starting atN L = exp 1 2αsb 0 exp [γ E ] that is present in the dQCD result. This contribution seems to be absent in eq. (97), but actually it is contained by the ratio of Γ-functions (eq. (95)). To see this, we set µ s = Q/N , as motivated by ref. [39,86]. We then find that where we have used The integral in eq. (104) can be recognized to give Substituting this result in eq. (104), we obtain With η a = A (1) a πb 0 ln(1 − 2λ) at the first non-trivial order, we indeed observe that this factor gives rise to the branch cut starting at λ = 1/2. Given that we can write where we set D Γ = 1 at LP LL and NLL, we have now established that where the first factor is contained in U (Q, µ s = Q/N ), while the second one is generated by the ratio of Γ-functions (eq. (95)). We can also relate D LP+NLP aa (which is equal to D LP aa up to NLP NNLL accuracy) to its SCET counterpart. From eq. (84), we see that this must be contained in −a γ S (Q 2 , µ 2 s ) + 2a γa (µ 2 s , Q 2 ). We first write Let us first examine the NNLL DY case, for which which may derived by adding the explicit forms of these coefficients as given in appendix D. We then find for eq. (110) We shall return to the role of this last term, but we first turn our attention to the Higgs production case, where the story is somewhat more involved. For the Higgs case we have We have an apparent discrepancy with the DY case, but must remember that U h (eq. (85)) has an additional prefactor that can be written as By substituting β 1 = (4π) 2 b 1 and β 0 = 4πb 0 , we observe that this factor cancels the first two terms of eq. (114), such that up to NNLL and for µ s = Q/N we have After these first two steps, we have established that up to NLL accuracy, and for µ s = Q/N By now truncating D Γ to LP NNLL instead (eq. (108) without the partial derivative towards N , which is of NLP), we obtain Note that D Γ gives an additional factor of ζ(2) at NNLL (see also eq. (43)), turning the 3ζ(2) discrepancy into a −ζ(2) difference between the wide-angle contribution D aa for dQCD, and the sum of anomalous dimensions for SCET. We will see that this remaining (last) term is resolved after including the hard function of SCET and dQCD. In eq. (118), we have achieved part of the LP result of eq. (97). What remains is to relate the hard function of dQCD, g 0 (α s ), to the product H(Q, Q)s a (0, Q/N ) in SCET (where we set µ s = Q/N , such that the second line of eq. (97) is 0). This comparison may be done easily. For DY we have (see eq. (157) of appendix D) for dQCD, and for SCET (second line of eq. (164)) For the Higgs case we have (eq. (158)) for dQCD, and for SCET (first line of eq. (164)) Note that in the SCET hard functions of eq. (120) and (122) there is an additional exponent with respect to the dQCD hard functions of eq. (119) and (121). These contributions precisely cancel the last term in eq. (118). By now combining these results in eq. (97), we find By comparing this result to the dQCD result of eq. (23) with h (1) set to zero, we see that the difference between dQCD and SCET at LP consists of the O(η a /N )-term. This term is included in the LP SCET result, but absent in the LP dQCD result. Therefore, we see that the difference between dQCD and SCET starts at NLP in Mellin space, as was noted before in ref. [39,81,85]. We now study how the NLP LL contribution modifies this result. First we see that η a (eq. (83)) at the first logarithmic order is given by Moreover, we set S a,LL (Q 2 /N 2 , Q 2 ) = − 2 αs g (1) (λ) using eq. (90). With this, we find that the NLP SCET result of eq. (98) can be written in Mellin space as where we have used the result of eq. (90). Adding this to the LP result restricted to LL (where H =s = 1, and a γ = 0), we find The NLP correction in the second factor is precisely the first-order expansion of the NLP exponent h (1) a (λ) of dQCD. Hence, at NLP LL, the contributions from SCET and dQCD are the same if one sets µ s = Q/N . Differences between the two formalisms originate from i) truncations of higherlogarithmic terms, ii) truncations on higher-power terms (i.e. a factor of 1/N 2 is included in dQCD once the NLP exponent is expanded up to O(α 2 s ) or beyond, while it is missing in the SCET formalism), and iii) running-α s effects.
We recall that we can obtain the NLP z-space contribution by deriving the LP LL result towards α s (µ s ) (eq. (92)). We now comment on how this compares to the dQCD result, where the NLP N -space result can be obtained via a derivative of the LP LL term towards N . The N -space form of eq. (92) is as these functions do not depend on z. Using the definition of the β-function and setting µ s = Q/N , the this equation can be rewritten as Since we have established that U aa,LL (Q, Q/N ) = exp 2 αs g (1) (λ) , we find that eq. (92) is equal to the derivative trick in dQCD (eq. (44)) at the first-order expansion of the exponent.
We have now established that the contributions from SCET and dQCD are the same at NLP LL accuracy if one sets µ s = Q/N . Two important remarks are in order. Firstly, for our SCET starting point we have included a factor of z −ηa in eq. (82). If we would not include this factor as in eq. (81), the ratio of Γ-functions in eq. (95) would instead be  The factor of η 2 a is removed by including z −ηa , which has a phase-space origin. In dQCD, we recognize this factor in the J k contribution in eq. (32), which also has a phase-space origin: it comes from the ln(z) expansion around z = 1 as performed in eq. (31). This factor is needed in dQCD to obtain the NLP contribution directly from performing a derivative on the LP contribution, as shown in eq. (36). However, even without performing the derivative trick, we truncate the dQCD NLP result at strictly NLP LL (resulting in the function h (1) a ). If one would omit the factor of z −ηa in the SCET resummation formula, one would introduce an NLP NLL term via the factor of η 2 a . Therefore, by including this factor, the SCET and dQCD formalisms are put on equal grounds. We shall see that the numerical impact of this factor is sizeable.
Secondly, although we have found analytical agreement at LP + NLP LL accuracy between dQCD and SCET , we will see that there are numerical differences, which especially grow sizable when the LP part is evaluated at NNLL. This is caused by the fact that in dQCD the NLP contribution is multiplied by the hard coefficient g 0 (α s ), whereas in SCET, the NLP LL result is added to the LP result. Indeed, in ∆ SCET,NLP , we do not see a hard function, so for SCET the NLP result will not get enhanced by constant contributions. This is simply the result of a choice, i.e. of including strictly NLP LL logarithms in z-space in the SCET NLP resummation formula eq. (89). As a remedy, we can upgrade the NLP SCET result for both DY and Higgs by usinḡ (130) We will explore the numerical consequences of this in the next section.

Numerical comparison of dQCD and SCET resummation at NLP
Here we show the numerical results of the dQCD and SCET comparison. For the LP NNLL + NLP LL dQCD results, we use expression eq. (49) in eq. (22), and we include g 0 (α s ) up to O(α s ). The SCET results are obtained by adding eq. (89) (the NLP result) to eq. (82) (the LP result with the factor of z −ηa ). We also show the results that are obtained excluding the factor of z −ηa from eq. (82), as in eq. (81), and those obtained by dressing the NLP SCET result by H(Q, µ h ) (eq. (130)). The results are shown in fig. 15. Interestingly, we observe that near-perfect agreement is found between dQCD and SCET at NLP LL accuracy if one includes both the hard function and z −ηa . If both these factors are included, we find that the difference between dQCD and SCET does not exceed In fig. 16 we examine how robust these results are under variation of the scale. We first examine the case where we vary µ between µ = 2Q and µ = Q/2 (top and middle panel). For the SCET LP result, we again use eq. (82) and for the NLP result, we use eq. (130). In general, the dQCD uncertainty is larger than the SCET scale uncertainty. Note that the dQCD scale uncertainty band does not contain the LP SCET central value for large values of Q. The situation is worse for the LP SCET scale uncertainty, where especially for the Higgs production case the LP dQCD central value is not contained in the uncertainty band. This changes at NLP. At NLP, the dQCD uncertainty band does not change in size, but with the SCET result lying much closer to the dQCD result, it is now contained inside the uncertainty band. For SCET, the uncertainty band does shrink at NLP, which is especially pronounced for Higgs production.
In dQCD, one implicitly varies the 'soft scale' via the inverse Mellin transform. If we equateN to Q/µ s , as we have done above, we readily see that all values of N (and hence of the soft scale) are accessible by integrating over N . In SCET one has to vary the soft scale explicitly, and we show the effect of this in the bottom panel of fig. 16. There, one may observe that the uncertainty obtained by varying the soft scale grows very large (especially for the Higgs case) once the NLP corrections are included. The soft-scale variation of the LP SCET result is small, and the LP dQCD result is not contained in this uncertainty band, however, the LP dQCD is contained in the soft-scale uncertainty band of the NLP SCET result.
In this section, we have compared the NLP LL SCET and dQCD contribution both analytically and numerically. At the analytical level, differences between NLP LL SCET and dQCD start at O(α 2 s ) and O(1/N 2 ). However, we find that the numerical impact of these terms is small, as we see near-perfect numerical agreement between dQCD and SCET for both the Higgs and DY cases. However, this agreement can only be obtained if one includes the hard function in the NLP SCET contribution at the same order as the LP SCET contribution, and if the factor of z −ηa is included in the LP SCET contribution. Without this latter factor, the difference between the NLP dQCD and SCET result starts at O(η 2 a /N ) = O(α 2 s ln 2N /N ), which is an NLP NLL contribution. This contribution turns out to be sizeable, as we have shown in fig. 15.

Conclusion
In this paper we have studied the role and impact of NLP corrections in colour-singlet production processes, with a particular focus on DY and single Higgs production. In section 2 we assessed the quality of the threshold expansion for Higgs and DY, for both the dominant (qq and gg) and the subdominant (qg) partonic channels. The threshold expansion of the dominant Higgs production channel is less well-behaved than that of DY, due to a substantial part of the Higgs partonic cross section arising from the small z region. The quality of the threshold expansion depends only marginally on the boson mass Q for DY, while the convergence for Higgs noticeably improves as Q increases. The threshold expansion of the off-diagonal qg-channel in Higgs production conver-gences only very slowly, whereas for DY convergence is already obtained after including the N 4 LP contribution.
We reviewed the resummation of leading-logarithmic NLP corrections in dQCD in section 3, and derived a slightly improved resummed expression involving a derivative with respect to the Mellin moment (eq. (45)). This we applied to the NLP LL + LP NNLL DY and Higgs cross sections (as well as di-Higgs and di-vector boson production, albeit at lower logarithmic accuracy) and performed numerical studies (section 4). For the DY/Higgs diagonal channels we found that the NLP LL resummation has a notable numerical impact on the NNLL resummed result matched to NNLO for both production processes, especially when compared to the effect of the N 3 LL contribution. A similarly sizeable impact was seen in the di-boson processes. The off-diagonal qg NLP LL contribution is more important for DY than for single-Higgs production given the large δ(1 − z) contributions in the latter process at NLO and NNLO, which enter the resummation formula via the matching coefficient g 0 and inflate the resummation effects in the Higgs diagonal channel.
In section 5, we performed an analytical and numerical comparison of the NLP LL resummed expressions in the dQCD and SCET frameworks. A derivative construction to obtain the NLP LL contribution, similar to the Mellin derivative in dQCD, can be employed in SCET (eq. (92)), where the derivative is performed in the coupling evaluated at the SCET soft scale µ s . At NLP we found from our analytical analysis that differences between the two formalisms originate from truncations of higher-logarithmic or higher-power terms, and running-α s effects. These differences resulted in very small numerical discrepancies once the NLP corrections were included; the two formalisms then agree numerically to O(0.5%)/O(1%) for Higgs/DY production after setting hard scales equal to Q, and the SCET soft scale equal to Q/N . To obtain this agreement, we argued that a factor of z −ηa is profitably included in the LP SCET expression, which originates from treating the kinematics of the scattering exact up to NLP. Indeed, sizeable numerical differences between the dQCD and SCET results of O(5 − 15%) in the DY case and O(15 − 40%) in the Higgs case are found if this factor is not included. Furthermore, we showed that the hard function of SCET that multiplies the NLP resummation should contain O(α s ) corrections as well, certainly for the Higgs case. Large numerical differences are found for the Higgs case without including this factor, where the δ(1 − z)-contribution, and hence the hard function itself, is large. For the DY case the need is somewhat less compelling, as its δ(1 − z)-contribution is smaller.
As part of our SCET-dQCD comparison we examined the factorisation/renormalisation-scale variation of the SCET and dQCD results, varying µ between Q/2 and 2Q, and the SCET soft scale µ s between Q/(2/N ) and 2Q/N . At LP, the dQCD and SCET central-value results do not lie in each other's scale-uncertainty bands. Conversely, at NLP, the central-value results show only a small numerical discrepancy and comfortably lie within each other's scale-uncertainty bands. The uncertainty bands do increase upon including NLP effects for both dQCD (by varying µ) and SCET (by varying µ s ). We expect this to be reversed once NLP NLL contributions would be included.
Our study of the numerical effects of NLP corrections in hadronic collisions is, we believe, a valuable addition to an area where most effort has so far been on the analytical side. It moreover validates these efforts by showing that NLP threshold corrections can have a notable impact, and should motivate further development of the understanding of NLP corrections.

Acknowledgments
EL and MvB acknowledge support from the NWO-I program 156, "Higgs as Probe and Portal", and MvB from the Science and Technology Facilities Council (grant number ST/T000864/1). JSD is supported by the D-ITP consortium, a program of NWO funded by the Netherlands Ministry of Education, Culture and Science (OCW). LV acknowledges support from the Fellini -Fellowship for Innovation at INFN, funded by the European Union's Horizon 2020 research programme under the Marie Skłodowska-Curie Cofund Action, grant agreement no. 754496. This paper is also based upon work from COST Action CA16201 PARTICLEFACE supported by COST (European Cooperation in Science and Technology).

A Fitted parton distribution functions
In this paper we rely on PDFs obtained from the LHAPDF library [90] and use the central member of the PDF4LHC15_NNLO_100 PDF set [54]. For our purposes these PDFs need to be converted to N -space to perform the resummation. To this end, we expand the PDFs on a basis of polynomials whose Mellin transforms may be computed directly. The functional form of the PDFs that we use is inspired by that used by the MMHT [91] collaboration, and reads xf (x) = A(1 − x) a 1 x a 2 (1 + by + c(2y 2 − 1)) + B(1 − x) a 3 x a 4 (1 + Cx a 5 ) for y = 1 − 2 √ x , (131) with 10 (real) fit parameters. We require that the fitted function lies within the 1σ error as given by the LHAPDF grid implementation of the PDFs in the entire domain. We have checked our set by comparing the fixed-order results obtained with the x-space form of the fitted PDFs with those obtained using the grid directly, and found that differences are smaller than the numerical integration error. A tabulated form in c++ format of the resulting fit parameters is available at [92].

B Normalization of the partonic cross section
In this appendix we discuss the normalization of the partonic cross section ∆ ij (z), with the particular aim to highlight the origin of the additional factor of 1/z appearing in Higgs production relative to DY. We start from the definition of the invariant mass distribution for DY, which according to the QCD factorisation theorem is given as whereσ ab (Q 2 /(x a x b S)) is the partonic cross section, and we drop the renormalisation/factorisation scale dependence for conciseness. For Higgs production an analogous equation holds, with the invariant mass distribution replaced by the total cross section. We first rewrite the partonic cross section in terms of z = Q 2 /(x a x b S) = Q 2 /s. This can be done by inserting where we recall that τ = Q 2 /S. With this, eq. (132) becomes This equation is then matched to eq. (14), which we repeat here The matching equation τ zσ ab (z) = σ 0 (Q 2 ) ∆ ab (z), is defined such that σ 0 (Q 2 ) contains the terms of the tree-level cross section which are z-independent, while ∆ ab (z) contains the z-dependent terms.
We are now in a position to perform the matching in eq. (136) for the DY and Higgs tree-level cross section. We start from the matrix element squared for the two processes, summed (averaged) over the final (initial) state partons .
The factor of s 2 in the Higgs matrix element squared is due to the derivative squared contained in the G a µν G a,µν H term of the effective Lagrangian. Given the one-particle phase space with the flux factor of 1/2s 1 2s dΦ 1 = d 4 q 2π δ (4) (p 1 + p 2 − q)δ + q 2 − Q 2 = π s 2 δ(1 − z) , the partonic cross section for the two processes readŝ We see that the DY cross section has an additional factor of z compared to Higgs production, whose origin is ultimately related to the dimensionful effective ggH vertex versus the elementary qqγ vertex. The coefficient ∆ At tree level the additional factor of 1/z is of course harmless, given the overall δ(1 − z). However, this factor is general, and it is present also at higher orders in perturbation theory, explaining the origin of the factor 1/z in eq. (9) (see also eq. (13)).

C Singular contributions at threshold (LP)
In order to asses how the parton flux weighs the partonic cross section, it is useful to consider the point-by-point multiplication of the partonic cross section and the parton luminosity factor L ij (τ /z)/z. For non-singular terms beyond LP this is trivial, but the singular contributions at LP consist of plus-distributions which have the non-local definition and T R = 1/2, C A = 3 and C F = 4 3 . The number of active flavours is denoted by n f and is set equal to 5 in this work. At O(α 2 s ), the solution to the β-function reads Alternatively, we occasionally use where β 0 = 4πb 0 , β 1 = (4π) 2 b 1 , etc.

D.1 dQCD
The initial state exponents for the LP LL (g (1) a ), NLL (g (2) a ) and NNLL resummations (g a ) are given by [62,100] g (1) a (λ) = where λ = b 0 α s lnN and α s ≡ α s (µ 2 R ). The N 3 LL function g (4) a (λ, Q 2 /µ 2 F , Q 2 /µ 2 R ) is extracted from the TROLL [40] code. The function h (1) which is added to account for the NLP LL terms is The coefficients A (n) a are given by [4,60,61] A (1) a = C a , with C q = C F and C g = C A . The perturbative (or hard) function g 0 (α s ) reads g 0 (α s ) = σ 0 1 + α s g 1 0 + α 2 While we extract the coefficients g i 0 from TROLL [40], we give their explicit form below for convenience. For DY we have σ 0 = σ DY , g DY,1 (3)  The other coefficients are given by