Fiducial Higgs and Drell-Yan distributions at N$^3$LL$^\prime$+NNLO with RadISH

We present state-of-the-art predictions for transverse observables relevant to colour-singlet production at the LHC, in particular the transverse momentum of the colour singlet in gluon-fusion Higgs production and in neutral Drell-Yan lepton-pair production, as well as the $\phi^*_\eta$ observable in Drell Yan. We perform a next-to-next-to-next-to-leading logarithmic (N$^3$LL) resummation of such observables in momentum space according to the RadISH formalism, consistently including in our prediction all constant terms of relative order $\alpha_s^3$ with respect to the Born, thereby achieving N$^3$LL$^\prime$ accuracy. The calculation is fully exclusive with respect to the Born kinematics, which allows the application of arbitrary fiducial selection cuts on the decay products of the colour singlet. We supplement our results with a transverse-recoil prescription, accounting for dominant classes of subleading-power corrections in a fiducial setup. The resummed predictions are matched with fixed-order differential spectra at next-to-next-to-leading order (NNLO) accuracy. A phenomenological comparison is carried out with 13 TeV LHC data relevant to the Higgs to di-photon channel, as well as to neutral Drell-Yan lepton-pair production. Overall, the inclusion of ${\cal O}(\alpha_s^3)$ constant terms, and to a lesser extent of transverse-recoil effects, proves beneficial for the comparison of theoretical predictions to data, leaving a residual theoretical uncertainty in the resummation region at the 2-5% level for Drell-Yan observables, and 5-7% in Higgs production.


Introduction
The experimental data collected in Run I and II at the Large Hadron Collider (LHC) has so far shown no significant deviation from the predictions of the Standard Model (SM) of particle physics. Since signals of new physics could emerge as tiny distortions in the spectra of sensitive observables with respect to the SM baseline, the availability of very accurate theoretical calculations, chiefly at the differential level, is of paramount importance. Processes featuring a colour-singlet system in the final state, such as Drell-Yan (DY) production or Higgs (H) gluon-fusion production, play a central role in the LHC precision programme. In particular, observables which depend only on the total transverse momentum of the associated QCD radiation represent an especially favourable environment both from the theoretical and the experimental viewpoint. On the one hand they feature comparatively low complexity, allowing one to push perturbation theory to its limits; on the other hand, their little sensitivity to multiparton interactions and non-perturbative modelling allows a particularly clean comparison between the theoretical predictions and the extremely precise experimental data, thereby challenging the accuracy of the former.
It is well know that fixed-order predictions must be supplemented with the all-order resummation of enhanced logarithmic contributions which arise in the phase-space region dominated by soft and/or collinear QCD radiation; by denoting with v a generic dimensionless transverse observable, i.e. one not depending on the radiation's rapidity (for instance p t /M or φ * η , M being the mass of the colour singlet), such a region corresponds to the v → 0 limit. The resummation of v spectra in colour-singlet production is customarily performed in impact-parameter b-space, where the phase-space constraints factorise [46,47]. Using the b-space formalism, the p t distribution in Higgs production has been resummed at next-to-next-to-leading logarithmic (NNLL) accuracy in Refs. [48][49][50], within the approach of [47,51], and in Ref. [52] using Soft-Collinear Effective Theory (SCET); N 3 LL resummation was considered in Ref. [53,54]. As for DY, p t and φ * η have been resummed in b-space at NNLL in Refs. [55][56][57][58][59] and at next-to-NNLL (N 3 LL) accuracy in Refs. [54,[60][61][62][63].
As an alternative to b-space resummation, the RadISH framework for the resummation of transverse observables in momentum space has been introduced in Refs. [64,65], which bases the resummation on a flexible Monte Carlo (MC) formulation (see also Ref. [66] for a study of directspace p t resummation in SCET). Resummed predictions at N 3 LL accuracy within the RadISH formalism have been presented for Higgs production at the inclusive level in Ref. [65] and within fiducial cuts in Ref. [67]. For Drell-Yan production, N 3 LL RadISH predictions for both p t and φ * η have been achieved in Refs. [67,68], and also considered in [69]. N 3 LL results for generic coloursinglet production, see for instance [70], are available through the automated MATRIX+RadISH interface [71,72]. Moreover, the momentum-space formulation is at the core of recent applications in the context of matching NNLO calculations with parton-shower simulations (NNLO+PS) [69,73,74].
In this article we consider again the Higgs p t distribution in gluon fusion, and the di-lepton p t and φ * η distributions in DY, and present state-of-the-art resummed predictions in which we consistently supplement known N 3 LL results with the inclusion of all constant terms of relative order α 3 s in the resummation, reaching so-called 'primed' accuracy N 3 LL . While the N 3 LO hard functions for DY and for Higgs production in the m top → ∞ limit have been known for some time [75][76][77], reaching N 3 LL accuracy for these processes requires, as also done in [12,34], to supplement the ingredients deduced in [78][79][80][81][82][83][84][85][86][87][88][89][90][91][92] with the quark and gluon transverse-momentum dependent (TMD) beam functions at N 3 LO, which were recently obtained via two independent calculations in Refs. [93][94][95]. Our predictions are further improved by the inclusion of transverserecoil effects, which we achieve by implementing in RadISH the prescription of Ref. [96].
We combine our resummed N 3 LL results with fixed-order differential spectra at NNLO accuracy from NNLOjet [14][15][16]38], and we present matched N 3 LL +NNLO predictions within fiducial cuts in comparison with 13 TeV LHC experimental data relevant to Drell-Yan di-lepton production [97], and to Higgs di-photon production [98].
This manuscript is structured as follows: in Sec. 2 we review the RadISH formalism for resummation in momentum space, up to N 3 LL order; Sec. 3 details the consistent inclusion of constant O(α 3 s ) terms, necessary to reach N 3 LL accuracy, and of transverse-recoil effects; in Sec. 4 we report on the tests we have performed to validate the correct implementation of the new contributions; phenomenological results at the LHC are presented in Sec. 5, and we give our conclusions in Sec. 6. We collect in Appendix A some formulae relevant for resummation up to N 3 LL , while Appendix B discusses subtleties related to the axial-vector structure of the three-loop DY form factor.

Momentum-space resummation in RadISH
The RadISH approach, developed in Refs. [64,65], is designed to resum recursively infrared and collinear (rIRC) safe observables [99] in momentum space. This is achieved by exploiting the factorisation properties of QCD squared matrix elements to devise a Monte Carlo formulation of the all-order calculation, effectively resumming large logarithms by generating soft and/or collinear radiation as an event generator of definite logarithmic accuracy.
The starting point is the cumulative probability for observable V ({p}, k 1 , . . . , k n ) (which, without loss of generality we assume as dimensionless) to be smaller than a certain value v where {p} =p 1 ,p 2 are the Born momenta of the incoming partons, and k 1 , . . . , k n are the momenta of radiated QCD partons. Even though the formalism is in principle extendible to generic rIRC safe observables, in the present article, as was done in Refs. [64,65], we focus on inclusive transverse observables: the former condition means V ({p}, k 1 , . . . , k n ) = V ({p}, k 1 + · · · + k n ), while the latter specifies that for a single soft emission k collinear to leg the observable can be parametrised as where M is the mass of the considered colour singlet, k t is the transverse momentum of k with respect to the beam axis, g (φ) is a generic function of the angle φ between k t and a reference direction n, orthogonal to the beam axis, d is a normalisation factor, and a > 0. For definiteness, the rescaled transverse momentum p t /M of the colour-singlet system features d = g (φ) = a = 1, while φ * η corresponds to d = a = 1, g (φ) = | sin(φ)|. In the soft limit, the cumulative cross section in (2.1) can be cast to all orders as where M is the renormalised matrix element for n real emissions (the case with n = 0 reduces to the Born contribution), [dk i ] denotes the phase space for the i-th emission with momentum k i , and the Θ function represents the measurement function for the observable under study. By Φ B we denote the Born phase space, while V(Φ B ) is the all-order virtual form factor relevant to the considered qq or gg reaction. The rIRC safety of the observable allows one to establish a well defined logarithmic counting for the squared amplitude [99,100], and to systematically identify the terms that contribute at a given logarithmic order. In particular, |M| 2 can be conveniently expanded in n-particle-correlated (nPC) blocks [65], defined as the contributions to the emission of n partons that cannot be factorised in terms of lower-multiplicity squared amplitudes. nPC blocks with higher n and loop order are logarithmically suppressed with respect to blocks with lower n and number of loops, so that an nPC block at l loops just enters at N n+l−1 LL accuracy.
The cumulative cross section in (2.3) contains exponentiated virtual IRC divergences in V(Φ B ), as well as real singularities in the multi-radiative squared matrix element. Such singularities are handled by introducing a resolution scale q 0 on the transverse momentum k t of radiation: rIRC safety ensures that blocks with total k t < q 0 , dubbed unresolved, contribute negligibly to the observable's value, and can be discarded in the evaluation of the measurement function; unresolved radiation thus exponentiates and regularises the divergences contained in V(Φ B ) at all orders. On the other hand, blocks harder than the resolution scale, referred to as resolved, must be generated exclusively, as they are constrained by the measurement function. The dependence of the prediction upon q 0 is guaranteed by rIRC safety to be power-like, hence the q 0 → 0 limit can be safely taken. For the observables considered in this paper, which solely depend on the total transverse momentum of QCD radiation, it is convenient to set the resolution scale to k t1 , where 0 < 1, while k t1 is the total transverse momentum of the hardest resolved block. We point out that the same resolution scale can be applied for the resummation of different observables, thereby allowing a flexible Monte Carlo implementation where multiple different resummations can be performed in a single framework, such as for instance the recently-introduced double-differential resummation of Higgs and leading-jet transverse momentum in gluon fusion [101]. After performing the above described set of operations, the all-order result for the cumulative cross section takes a particularly compact form in Mellin space, where convolutions with parton densities reduce to algebraic products. We introduce Mellin moments of generic functions g(z) as , and define f as the array containing the 2n f + 1 parton densities, (n f being the number of light flavours), whose DGLAP [102][103][104] evolution between scales µ 0 and µ reads with P the path-ordering symbol,P f (a)f (b) the regularised collinear splitting functions, and f (a) the flavour of the a-th entry of f . For notational simplicity, for the time being we consider only flavour-conserving kernels, so to make the Γ matrix diagonal and drop the path ordering; we will relax this assumption by the end of the section. The cumulative cross section differential in the Born variables can be written, with the convention of [65], as where |M B | is the Born squared matrix element, the sum runs over all allowed Born flavour combinations, Ω B denotes a set of internal phase-space variables of the colour-singlet system, and the integration contours C 1 and C 2 in the double inverse Mellin transform lie along the imaginary axis to the right of all singularities of the integrand. TheΣ matrix encodes the effect of parton-density DGLAP evolution from scale µ 0 , as well as that of flavour-conserving radiation evolving the partonic cross section. For inclusive observables, its all-order expression under the above assumption on Γ N is 1 The last two lines of eq. (2.6) reduce to Θ (v − V ({p}, k 1 )) for n = 0.
H(µ R ) represents the finite contribution to the virtual form factor, evaluated at the renormalisation scale µ R , and has a perturbative expansion of the form It satisfies a flavour-conserving renormalisation-group evolution equation stemming from the running of its coupling where we unambiguously dropped the c index in C (n) and in Γ (C) , for the sake of brevity. The R function encodes the contribution from radiation off leg which conserves the momentum fraction of the incoming partons and the flavour c of the emitter, namely [R ] ab = R δ ab . It is related to the Sudakov radiator R, with entries [R] ab = R δ ab , by Finally, the anomalous dimensions A and B encode the inclusive probability |M(k)| 2 inc [65] for a correlated block of arbitrary multiplicity to have total transverse momentum k t ; they admit a perturbative expansion as The structure of (2.6) shows the different contributions of resolved and unresolved radiation. The former, encoded in the third to fifth line, is represented by an ensemble of emissions (more appropriately: of correlated blocks treated inclusively) harder than k t1 , with contributions from flavour-diagonal radiation as well as from exclusive DGLAP-evolution steps. Conversely, the exponentiated unresolved emissions combine with the all-order virtual form factor giving rise to the Sudakov exponential e −R( kt1) . The factor C c1;T N1 (α s (µ 0 )) H(µ R ) C c2 N2 (α s (µ 0 )) encodes the hardvirtual corrections to the form factor, and the collinear coefficient functions. The coupling of the latter is evaluated at scale µ 0 and subsequently evolved inclusively up to k t1 by the operator containing Γ (C) N in the second line of (2.6). Similarly, the parton densities are DGLAP-evolved from µ 0 up to k t1 by Γ N .
As shown in Ref. [105], for gluon-fusion processes the structure in eq. (2.6) must be supplemented with the contribution from the (flavour-diagonal) G collinear coefficient functions, describing the azimuthal correlations with initial-state gluons. This contribution, starting at O(α 2 s ), i.e. N 3 LL order, is included in the above formulation by adding to eq. (2.6) an analogous term where one performs the replacements In the following, this contribution is understood whenever not explicitly reported.
The evaluation of eq. (2.6) at this point may be simplified by exploiting again rIRC safety. The latter ensures that the transverse momenta of all blocks in the resolved ensemble are parametrically of the same order, as blocks that are significantly softer than k t1 do not contribute to the evaluation of the observable and are accounted for in the radiator. All resolved contributions in eq. (2.6) with argument k ti ≡ ζ i k t1 can thus be Taylor-expanded about k t1 , with subsequent terms in the expansion being more and more logarithmically suppressed, since ζ i is of O(1). Analogously, unresolved quantities depending on k t1 can be expanded about k t1 : the ensuing logarithms ln(1/ ) exactly cancel the logarithmic -dependence of the corresponding terms in the resolved radiation, conveniently achieving an all-order subtraction of IRC divergences.
Aiming for N 3 LL accuracy, one needs to retain only the following terms in the Taylor expansion of the unresolved quantities as well as of the resolved contributions, which are suppressed by one logarithmic order with respect to the corresponding unresolved ones: where R (j) (k t1 ) = d j R(k t1 )/dL j , L = ln(M/k t1 ), and the ellipses denote neglected N 4 LL terms.
The loop expansion of the involved anomalous dimensions obeys an analogous perturbative counting. A further significant simplification stems from the fact that, at a given logarithmic accuracy, one needs to retain subleading terms in the above expansions only for a limited number of resolved blocks. For instance, at N k LL, only up to k − 1 resolved blocks need to feature a ln 1/ζ i correction in R , as the simultaneous correction of k factors of R affects N k+1 LL. Unresolved contributions are expanded correspondingly, in order to cancel the divergences of the modified resolved blocks to the given logarithmic order. By means of the above expansions, the master formula (2.6), which is valid to all logarithmic orders, at N 3 LL (and, as we will show in the next section, at N 3 LL as well) reduces tô (2.14) The final operation is to rewrite eq. (2.14) in direct (as opposed to Mellin) space, which requires little effort at this point, as a very limited number of exclusive evolutions steps have been retained in the above expression. In particular, at N 3 LL, only up to two hard-collinear resolved emissions are needed, and one can relax the above assumption of flavour-conserving real radiation by including flavour-changing kernels in the DGLAP-evolution contributions in momentum space. This amounts to the following identifications, valid at N 3 LL: where we defined ∂ L ≡ d/dL, L = ln(M/k t1 ), β 0 is the lowest-order contribution to the QCD beta function, L is the parton luminosity (see Appendix A for its explicit expression at the various logarithmic orders, and Sec. 3.1 for a discussion about standard and improved luminosities in the context of N 3 LL -accurate predictions), and After expressing the logarithms of 1/ as dummy radiative integrals [100] according to and introducing the average of a function G({p}, {k i }) over the measure dZ in which the dependence upon the regulator exactly cancels to all orders, one finally gets at N 3 LL where the explicit factors of α s (k t1 ) are defined as α s (k t1 ) = α s /(1 − 2α s β 0 L), and α s = α s (µ R ) unless stated otherwise.
We conclude this review of the RadISH approach with two remarks on the master formula (2.18). First, we note that the logarithms resummed there are of the form L = ln(M/k t1 ). It is convenient to introduce the resummation scale Q, of order M , as an auxiliary scale to be varied in order to probe the size of neglected logarithmic corrections, and resum logarithms ln(Q/k t1 ). This is formally achieved by splitting L = ln(Q/k t1 ) + ln(M/Q), by assuming the hierarchy ln(Q/k t1 ) ln(M/Q), valid in the IRC limit, and by expanding L around ln(Q/k t1 ) at the relevant logarithmic accuracy. Second, when the resummed results are matched to a fixed-order prediction, it is desirable to enforce the former to vanish in the hard region k t1 Q of the v spectrum, reliably described by the latter. This can be achieved by modifying the resummed logarithms ln(Q/k t1 ) by means of power-suppressed terms, negligible at small k t1 . A possible choice for modified logarithmsL is where p is a positive real parameter chosen so that the resummed differential spectrum vanishes faster than the fixed-order one at large v. The above prescription induces a jacobian J (k t1 ), which ensures the absence of subleading-power corrections with fractional α s powers in the final distribution, still keeping the k t1 → 0 region unmodified. We stress that the procedure of logarithmic modification is not just a change of variables, as it does not affect the observable's measurement function. As a consequence, the final resummed result shows an explicit p dependence through power-suppressed terms, which however, after matching, will cancel up to the accuracy of the fixed-order component.
In the following developments of the article we understand the procedure of logarithmic modification, which formally corresponds to considering the logarithmic region k t1 < Q and working in the p → ∞ limit of (2.19) and (2.20); moreover we redefine L ≡ ln(Q/k t1 ), in order not to unnecessarily clutter our formulae.

Consistent inclusion of N LL and recoil effects
In this section we discuss how the formalism detailed above can be upgraded to N 3 LL accuracy, which amounts to supplementing the N 3 LL result with the complete set of constant contributions of relative order O(α 3 s ) with respect to the Born. Such contributions formally pertain to the logarithmic tower α n s L n−3 , namely they are a subset of the N 4 LL correction, however they are of particular relevance since their inclusion suffices for the perturbative expansion of the resummed cumulative cross section Σ(v) to correctly encode all terms of order α n s ln 2n−6 (1/v). The definition of 'primed' accuracy requires to specify more precisely how these constant terms are actually included in Σ(v). In particular, at N k LL order, all choices leading to differences beyond α k s and N k LL accuracy are legitimate, such as, for instance, the argument of the coupling constant multiplying the highermost-order coefficient functions present in the 'primed' luminosity factors. In this work we adopt as our default 'primed' predictions those obtained by evaluating such a coupling at the scale k t1 of the hardest emission, which is the correct scale one would have to use for N k+1 LL accuracy. In Sec. 3.1 we will discuss this choice in more detail, and assess its impact in Sec. 5. We stress that, although in other formalisms 'primed' accuracy may correspond to encoding different subleading contributions with respect to ours, we decided not to introduce a new nomenclature for our improved predictions, given that, regardless of the formalism, N k LL results are anyway designed to upgrade N k LL ones by the inclusion of the whole set of constant α k s terms. The resummation formula presented in eq. (2.6) is formally valid to all logarithmic orders. However, the accuracy of its practical implementation is limited by the fact that the quantities it features, such as anomalous dimensions and coefficient functions, are known to finite perturbative order, and by the fact that, for computational convenience, the expansions detailed in eq. (2.12) and (2.13) have been performed to arrive at expression (2.14) in Mellin space and (2.18) in momentum space. Achieving full N 3 LL amounts to lifting the subset of such approximations that affect thirdorder constant contributions.
Focusing on the structure of eq. (2.14), and recalling that the weight of the hardest resolved radiation k t1 provides at least one power of α s , it is immediate to verify that the inclusion of further logarithmic derivatives in the exponent of the second line, as well as in the resolved ensemble, only affects O(α 4 s ) terms. We conclude that the structure of eq. (2.14) is sufficient as is to achieve N 3 LL accuracy: one just needs to evaluate its contributions to appropriate perturbative order, and to upgrade the conversions (2.15) to momentum space, which we address in turn in the next subsections.

O(α 3
s ) constants from radiator, hard, and coefficient functions The first source of constant O(α 3 s ) terms we consider is the radiator defined in eq. (2.9), which can be rewritten as: where L = ln(Q/k t1 ), L M = ln(Q/M ), λ = α s β 0 L, λ M = α s β 0 L M , and α s = α s (µ R ). The g k functions encode the resummation of N k−1 LL logarithmic towers α n s L n+2−k ; they are explicitly reported in Appendix B of [65] for k ≤ 4. From the above integral expression one notes that all constant O(λ 0 ) terms in the radiator vanish as λ M = 0, i.e. they are proportional to powers of L M . By introducing the expansion of the g k functions in powers of λ as g k (λ) = ∞ n=0 g k,n λ n , with g 1,0 = g 2,0 = 0, the constant g k,0 can be inferred by solely analysing g j functions with j < k: this stems from the fact that g k,0 is responsible for the cancellation of a well-defined part of the Q-dependence in the N k−2 LL radiator. At LL one has where the Q-dependence starts at NLL order, in the coefficient of the α n s ln n (M/k t1 ) term. The latter dependence is compensated by including g 2 in the radiator, so that the L g 1 + g 2 sum features a Q-dependence starting at NNLL, in the coefficient of α n s ln n−1 (M/k t1 ). In turn, the Q-dependence in the α s ln 0 (M/k t1 ) term of the L g 1 + g 2 sum is exactly compensated by the g 3,0 constant: from which g 3,0 can be read off. We note that there is no other constant term contributing to g 3,0 since, as observed above, all constants g k,0 are proportional to non-zero powers of L M . By generalising this argument, the analysis of the α k−2 s ln 0 (M/k t1 ) term in the radiator including up to g k−1 allows to deduce g k,0 , yielding the all-order expression 2 In our approach the constant part of the radiator, containing up to g 5,0 in the N 3 LL case, is then expanded in powers of α s , which avoids the presence of any exponentiated constants. Such terms are included in the hard-virtual function H contributing to the luminosity at the various perturbative orders, which thereby acquires an explicit Q-dependence, as detailed in Sec. 3.2.
Further O(α 3 s ) constant terms in the one-emission contribution to eq. (2.14) originate from the use in RadISH of anomalous dimensions A CSS, , B CSS, calculated for a b-space resummation, where the Sudakov radiator is usually defined as The conversion between the Heaviside and the Bessel function is absorbed into a redefinition of A CSS, , B CSS, , H CSS , and C CSS by means of the relation [106] 1 − J 0 (bk t ) = 1 + ζ 3 12 which starts being non-trivial at N 3 LL [65], involving the third logarithmic derivative R of the LL radiator function g 1 . In order to incorporate O(α 3 s ) constant effects, it is necessary to extend this construction to the third derivative of function g 2 , as well as to include the interference of the one-loop hard function H (1) with the third derivative of g 1 . This results in a constant term is the strong-coupling order of the Born squared amplitude (e.g. d B = 2 for Higgs production, and d B = 0 for Drell-Yan production). We include the δH (3) constant in the three-loop hard-virtual function H (3) , see Sec. 3.2. The conversion described above for the Sudakov radiator applies analogously to the partondensity and to the coefficient-function evolution exponent in eq. (2.6). While the third derivatve of the latter starts contributing at O(α 4 s ), the former generates an O(α 3 that we include into the third-order coefficient function C ij , see Sec. 3.2. Finally, we note that subleading terms in eq. (3.7) are proportional to the fifth logarithmic derivative of the Sudakov radiator, hence they start contributing at O(α 4 s ), and are neglected in this article.
) of eq. (2.6) are another source of constant terms, included in the luminosity factors of eq. (2.18). The latter admit a perturbative expansion that, in turn, originates from the ones of the hard-virtual function H, and of the collinear coefficient functions C and G. Such expansions, already introduced in Sec. 2, are reported here for convenience, using explicit flavour indices: where µ is the same scale at which parton densities are evaluated, and µ R is the renormalisation scale. In eq. (3.10) we only retained the perturbative orders needed to assemble an N 3 LL -accurate luminosity, where for the first time one needs the third-order coefficient and hard functions C (3) and H (3) , respectively, and the second-order azimuthal coefficient function G (2) , as discussed in the following. At N k LL order, the luminosity in the first line of the RadISH master formula (e.g. L N 3 LL (k t1 ) in eq. (2.18)) contains all constant O(α k−1 s ) terms, which properly pertain to the N k LL logarithmic tower, whence the subscript labelling L. The coupling constants of the involved C and G coefficient functions have to be evaluated at the same scale at which the parton densities are evaluated, i.e. µ = k t1 . On the other hand, when working at N k LL accuracy, one has the freedom to choose whether the O(α k s ) constant terms included in L N k LL (k t1 ) are evaluated with a fixed scale, e.g. µ = µ R , or a running scale, for instance µ = k t1 : this ambiguity only affects terms starting from O(α k+1 s L), namely non-constant N k+1 LL contributions beyond accuracy.
The above discussion can be easily illustrated focusing on the lowest order at which it applies, namely NLL : the NNLL luminosity reads whereas, at NLL , one is allowed to define either L NLL = L NNLL or Other choices for the running coupling of the O(α k s ) terms are of course equally allowed, and we consider these two as representative of the genuine perturbative ambiguity underlying 'primed' predictions. We refer to results obtained with these two different choices as with or without running coupling, respectively. At any order above NLL , L N k LL (k t1 ) features a similar ambiguity in constant terms of order O(α j s ), with j < k, where one is allowed to run the coupling at arbitrary loop order, provided the latter is ≥ k −j. For consistency, in the L N k LL (k t1 ) luminosity without running coupling, constant O(α j s ) terms are evolved at k − j loops, while with running coupling they are evolved at k − j + 1 loops. Similar considerations apply to the luminosity factors appearing in the contributions with one and two special emissions, i.e. the lines beyond the first in eq. (2.18).
The full expressions for the upgraded luminosities up to N 3 LL , with and without running coupling, are given in Appendix A. For the phenomenological N 3 LL presented in this paper, we have considered both options, choosing as our default the one with running coupling, as it includes a whole tower of correct N k+1 LL effects. We will show the effect of this choice quantitatively in Sec. 5.

Extraction of hard and collinear coefficient functions at O(α 3
s ) In this subsection we discuss the extraction of the hard-virtual function H, and of the collinear coefficient functions C and G needed for N 3 LL accuracy.
The H (n) coefficient of the hard-virtual function in eq. (3.10) is obtained from the quark and gluon form factors at n loops. Except for ζ 3 -contributions analogous to (3.8), it coincides with the nth term of the perturbative expansion of is the Wilson coefficient of the ggH effective vertex in the MS scheme [75,76] and C(α s (M ), M 2 , M 2 ) is the hard matching coefficient of Ref. [77], evaluated with time-like virtuality. Explicit expressions for n = 1, 2 are reported in eqs. (3.30) and (3.31) of Ref. [65]. The third-order hard coefficients for Higgs and Drell-Yan production read where M is the mass of the colour singlet, and L t = ln(M 2 H /m 2 top ), being M H the Higgs mass. The above expression for H (3) q (M ) matches eq. (7.8) of Ref. [77], without the term proportional to N F,V : this term originates from the structure of the vector and axial-vector couplings of the neutral Drell-Yan process, and its presence implies that the Born matrix element cannot be exactly factored out of the hard-virtual coefficient at two and three loops. We discuss the physical reasons for this subtlety, and how we handle it, in Appendix B.
As discussed in Sec. 3.1, in the RadISH formalism the shift δH (3) defined in eq. (3.8), as well as the constant parts of the radiator in eq. (3.4), are absorbed in the hard coefficient H (3) . After taking into account the explicit dependence of the hard function upon the renormalisation scale, the final expression for H (3) Obviously, in eq. (3.14), the dependence of the hard coefficient H (n) (M ) on the process at hand is understood, as is the case for the anomalous dimensions contained in g k,0 . The above expression matches exactly what we implemented in the RadISH code.
The collinear coefficient functions C (3) are extracted from the transverse-momentum-dependent (TMD) parton-density functions (PDFs). These, in turn, are obtained combining the TMD beam functions and the soft function, computed at third order in [93][94][95] and [83,84], respectively. The G (2) function can be instead extracted from the computation of the linearly-polarised gluon TMD PDFs at two loops [88,107]. Since our starting points for the extraction of C (3) and G (2) are Refs. [83,88,93,95], in order to make contact with the notation used therein, we recall that, when expressed in terms of TMD beam and soft functions, the factorisation formula for transversemomentum resummation in b-space has the schematic structure 3 where H, B, and S ⊥ are the hard, beam, and soft functions, respectively. Although the exact correspondence between the RadISH formalism and resummation in impact-parameter space has been discussed elsewhere [65], by comparing eq. (2.6) with eq. (3.18) it is easy to see that the C and G functions are to be extracted from the combination As explained in Refs. [88,93,95,108] (for instance, eq. (3.8) of Ref. [95]) this is also the combination that allows one to define a TMD parton density that does not depend on either the rapidity regulator used for the computation of B and S ⊥ , or the rapidity scale.
The beam function B admits an Operator Product Expansion (OPE) onto the collinear PDFs. After renormalisation, and after the remaining collinear divergences are reabsorbed into the collinear PDFs, the unpolarised quark and gluon beam function are defined by the coefficients of the OPE, i.e. by the so called perturbative matching coefficients [93] for the quark case, and in eq. (3.4) of Ref. [95] for both the quark and the gluon cases. In addition to the I gi (ξ, b ⊥ ) coefficient, the tensor structure of the gluon beam function contains a second term, associated to a linearly polarised gluon, whose coefficient is denoted by I gi (ξ, b ⊥ ) (see eq. (2.4) of [95]): such coefficient gives rise to the G gi collinear function. 4 The results of Refs. [88,93,95] and of Ref. [83] contain the boundary conditions for the TMD soft and beam functions, as well as the complete scale-dependent expressions for B and S ⊥ , obtained solving their evolution equations. We are not interested in the latter, as the evolution of the coefficient functions is obtained in RadISH by means of the anomalous dimensions Γ (C) and Γ (G) . According to the above discussion, the C where the I  [83]), whose overall colour factor is C F (C A ) for Drell-Yan (Higgs) production, are the boundary conditions for S ⊥ . Similarly, the G  [88,95], and we have inserted them in the RadISH code, following the conventions of Refs. [110,111] for the flavour-decomposition of coefficient functions. As a cross check, we have verified that we obtain for C (2) ij (z) the same result we extracted in Ref. [65]. The C (3) ij (z) expressions contain harmonic polylogarithms (HPLs) of weight up to 5, which we efficiently evaluate via the fortran routine hplog5 [112]. 5 As far as the numerical implementation is concerned, we perfectly reproduced Fig. (2) of [93], and we verified that our fortran implementation matches the numerical results obtained using Mathematica and the package HPL [113]. Moreover, we also checked our implementation by comparing against the N 3 LO TMD PDFs results obtained in Refs. [94], finding perfect agreement.
As discussed in Sec. 3.1, in the RadISH formalism we absorb in the third-order coefficient function the shift δC (3) ij (z) defined in eq. (3.9). Furthermore, in order to match our resummed results to fixed-order calculations that feature α s (µ R ) and f (µ F ), we write the factors of α s (k t1 ) and f (k t1 ) appearing in the luminosities in terms of α s (µ R e −L ) and f (µ F e −L ), respectively, with L = ln(Q/k t1 ), absorbing the ensuing constant difference in the coefficient functions. This gives rise to an explicit µ F,R dependence in the latter, which is reported in eq. (4.6) of Ref. [65] for C (1) ij and C (2) ij . As for C (3) ij , we document such a dependence in the following equation: Analogously, for G gj (z) one has The above equations exactly match the expressions implemented in the RadISH code. In order for the next section to be notationally consistent with the previous ones, we will still denote parton densities and coupling constant as f (k t1 ) and α s (k t1 ) in the following formulae, understanding the above discussion.

O(α 3 s ) constants from multiple resolved emissions
We now turn to the description of O(α 3 s ) terms stemming from the two-emission contribution to (2.14). The first immediate correction comes from the last identification in the list (2.15), where we now need to retain Γ where we defined The following correction features for the first time the contribution of Γ (C) : where we have used the evolution equation (2.8) to evaluate Γ (C,1) = −4πβ 0 C (1) , and Next, an O(α 3 s ) contribution coming from the derivative of the DGLAP anomalous dimension in the last line of (2.14) reads Analogously, a constant O(α 3 s ) term is induced by a luminosity upgrade L NLL (k t1 ) → L NLL (k t1 ) in the fifth line of eq. (2.18), where L NLL (k t1 ) was introduced in Sec. 3.1. Conversely, the three contributions are already accounted for by the third line of eq. (2.18), and need not be added.
The final terms to be considered are corrections to the three-emission contributions. They feature a term with three lowest-order DGLAP-evolution matrices and two terms with the second derivative of the radiator and Collecting all contributions, our final formula for direct-space resummation at N 3 LL reads We stress that the comments on the modified logarithms and jacobian factor reported below eq. (2.18) apply unchanged to eq. (3.33) as well.

Transverse-recoil effects
In order to realistically simulate the kinematics of the singlet's decay products, we have implemented in our framework the default transverse-recoil prescription of [96] to account for the singlet recoiling against initial-state QCD radiation. The procedure amounts to considering the differential spectrum with respect to observable v, and to boosting its underlying Born kinematics from a rest frame of the singlet (specifically, the Collins-Soper one [114] in the default prescription) to the laboratory frame: there the singlet has transverse momentum equal to q t (v), where q t (v) = M v (or q t (v) = M v/| sin φ|, with φ the singlet's azimuthal angle) if v = p t /M (or v = φ * η ). Fiducial selection cuts are then applied on the boosted Born kinematics.
As argued in [63], see also [54], the inclusion of recoil effects via the prescriptions of [96] is sufficient to account for all linear power corrections in presence of fiducial cuts, together with their resummation with the same accuracy as the leading-power terms, for observables which are azimuthally symmetric at leading power, such as p t /M . Let us briefly discuss the technical implementation of recoil effects in the RadISH code. For each m-parton contribution to eq. (3.33), as defined by the Θ(v − V ({p}, k 1 , . . . , k m )) measurement functions, we evaluate the transverse momentum q t (v) of the colour singlet and its azimuthal angle φ, and we apply the above mentioned boost. In order to enforce fiducial cuts on the boosted Born system, we modify each measurement function in eq. (3.33) as where the dependence on k 1 , . . . , k m in Θ cuts encodes the effect of the boost (i.e. Θ cuts equals 1 or 0 if the boosted Born configuration passes or not the cuts). On the contrary, in absence of recoil effects, the action of the cuts does not depend on momenta k 1 , . . . , k m : the constraint Θ cuts (Φ B , {k 1 , . . . , k m }) reduces to Θ cuts (Φ B ) and factorises out of the resummation formula, therefore eq. (3.33) is calculated only for the points which pass the fiducial cuts. Finally, in order to match the resummed result with fixed-order predictions, when transverserecoil effects are included we also need to modify the perturbative expansion of the resummation. As detailed in Ref. [65] (see in particular Sec. 4.2), in the default code the latter expansion is computed at the cumulative level, and expressed as a combination of classes of 'master' integrals. Since the recoil procedure entails boosts on the differential spectrum, we now first compute the derivative of the expansion at a given value v, and then apply fiducial cuts on the boosted kinematics, consistently with what is done in the resummation component.

Validation
In this section we discuss the tests we performed to validate our implementation of N 3 LL effects in the RadISH code.
A first robust check is achieved by comparing the α 3 s expansion of the momentum-space resummation formula for p t with the analogous expression derived starting from the cumulative p t cross section in b-space: where J 1 is the second Bessel function, and R b is the radiator as written in (3.5) in terms of anomalous dimensions A CSS, , B CSS, . We stress that this test has the virtue of allowing to assess at the analytic level the correctness of the δH (3) and δC (3) ij terms derived in Sec. 3.1. The inverse Fourier transform (4.1) can be calculated by Taylor-expanding the radiator and the luminosity factor around b = b 0 /p t at the appropriate order. This allows to write the cumulative cross section as where R CSS was introduced in eq. (3.5), and c n (p t ) are coefficients encoding luminosity and radiator information. The integrals in eq. (4.2) are then readily obtained as derivatives with respect to R CSS (p t ) of the generating functional This procedure provides an analytic expression to be directly compared with the momentumspace expansion, which is written (see the discussion in Sec. 4.2 of Ref. [65]) as a linear combination of classes of 'master' integrals. Such master integrals are evaluated analytically up to O(α 2 s ) 6 , while we resorted to high-accuracy numerical integration for those entering at O(α 3 s ). By comparing the two expressions for each relevant combination of A (n) and B (n) anomalous dimensions, and retaining full renormalisation, factorisation, and resummation scale dependence, we achieved complete analytic agreement at order α 2 s , and numerical agreement at or below the permyriad level for all terms entering at α 3 s , which is the numerical accuracy level of the master integrals. An analogous check was performed in the case of the φ * η expansion, finding similar agreement. As a further stringent test, we have numerically checked that the µ R , µ F , and Q dependence of our N k LL cumulative results cancels exactly at order O(α k s ), and is of relative order O(α k+1 s L) with respect to the Born, i.e. a pure N k+1 LL effect. In order to perform this test, we evaluate our expressions in the small-coupling regime α s 1 with a set of analytic toy PDFs [115]. If the dependence on e.g. the renormalisation scale µ R is implemented correctly, one must obtain where d B = 2 (d B = 0) for Higgs (DY) production, and λ is an O(1) rescaling factor. For sufficiently small α s values, we have tested the exact cancellation of the scale dependence by confronting the ratio ∆ N k LL (v; κα s )/∆ N k LL (v; α s ) against its expected scaling κ d B +k+1 . A similar test has been successfully performed on the expansion of the resummation formula in powers of α s . We have also explicitly checked that the artificial introduction of small bugs in the coefficients of the scaledependent terms results in clearly visible violations of the test, whose successful outcome then strongly corroborates the robustness of our implementation. Finally, as an internal self-consistency test, we compare the resummed result for Σ N k LL (v) to its O(α k s ) expansion in the asymptotic v 1 limit. Owing to the presence of modified logarithms, the two expressions are expected and numerically checked to coincide in such a limit, which also ensures the absence of residual exponentiated constants in the resummed expressions.

Phenomenological results at the LHC
In this section we present predictions up to N 3 LL +NNLO 7 relevant for neutral Drell-Yan leptonpair production, and for gluon-fusion Higgs production and decay to a photon pair, at the 13 TeV LHC. For both processes we consider inclusive and fiducial setups, the latter allowing a direct comparison with experimental data, without relying on Monte Carlo modelling for acceptances. We stress that the availability of theoretical results at the fiducial level is guaranteed by the fact that our resummmation formalism is fully differential with respect to the Born phase-space variables.
In principle, the availability of an N 3 LL resummation would allow us to obtain results for the N 3 LO fiducial Drell-Yan and Higgs cross sections by means of a slicing technique such as q Tsubtraction [116]. It is however well-known that, especially in presence of symmetric cuts on the p t of the singlet's decay products, such a technique requires to push the slicing parameter down to very small values, requiring an extreme control on the stability of the numerical calculation in the far IRC regime. This in turn translates into the necessity of dedicated high-statistics fixed-order predictions, to minimise possible numerical fluctuations. We thus refrain from quoting fiducial cross sections at N 3 LO in this article, and leave this development for future studies.
Aiming at reliable predictions across the entire v phase space, we match our resummed results with fixed-order differential spectra computed with the NNLOjet code, and used in previous works [67,68]. The matching is designed to reproduce the resummed prediction in the v → 0 region, dominated by soft/collinear QCD radiation, while reducing to the fixed-order calculation in the hard tails v 1. In [67] we adopted a multiplicative matching at the cumulative level. Besides an improved numerical stability in the v → 0 limit, where the cancellation between the fixed-order result and the perturbative expansion of the resummation can be delicate, a cumulative multiplicative scheme had the advantage, for processes with known total N 3 LO cross section, of extracting the constant O(α 3 s ) terms from the fixed-order result trough matching. Since such terms are now included directly in the resummation at N 3 LL accuracy, a multiplicative scheme is no longer advantageous in this particular respect. Moreover, as discussed in Sec. 3.4, we now include in our framework a transverse-recoil prescription to improve the kinematical description of the singlet's decay products, which is implemented at the level of the differential v spectrum.
For these reasons, in this phenomenological study we adopt as our default a differential matching belonging to the additive family, defined as /dv represents the perturbative expansion of the resummed spectrum at the same order. The Z(v) factor, that we choose as [65] is designed to enforce a dampening of the resummation component in the hard region of the spectrum, while leaving the v → 0 limit unaffected. We set u = 2 (we stress that u must be > 1 not to induce linear power corrections), and h = 3 as our defaults; we take a central v 0 = 1 , and consider a variation of v 0 in the range [2/3, 3/2] around its central value in order to reliably estimate matching systematics. Our reference value for the parameter p appearing in the definition (2.19) of modified logarithms is p = 4; we have checked that a variation of p by one unit does not induce significant differences. We also present results obtained through a multiplicative matching at the differential level, defined as where Z(v) is the same function introduced for the additive matching. Analogously to the additive case, the matching in eq. (5.3) only acts at the level of quadratic power corrections for u = 2.
We stress that all of our resummed calculations feature a Landau singularity arising from configurations where QCD radiation takes place at transverse-momentum scales k t ∼ M e −1/(2β0αs(M )) ∼ 0.1 GeV. In the predictions we present in the following, we set our results to zero when the hardest radiation's transverse momentum is below the singularity. This prescription has a negligible impact on differential spectra for typical values of M as, due to the vectorial nature of the considered observables [46,65], the v → 0 limit is dominated by radiation at the few-GeV scale, significantly harder than the Landau scale. We however stress that for a precise description of this kinematic regime, a thorough study of the impact of non-perturbative corrections, not included in the present article, would be necessary.
We finally recall that in all predictions shown in the following we adopt the NNLO DGLAP evolution for parton densities. Although the NNLO corrections to the evolution are formally of N 3 LL order, we include them also in the NLL and NNLL results to ensure an identical treatment of heavy-quark thresholds. Parton densities are evolved from a scale µ 0 ∼ 1 GeV upwards by means of the Hoppet package [117], which is used as well to handle all parton-density and coefficient-function convolutions.

Drell-Yan results
For Drell-Yan phenomenology, we consider pp collisions at 13 TeV centre-of-mass energy, and we use the NNLO NNPDF3.1 PDF set [118] with α s (M Z ) = 0.118 through the LHAPDF interface [119]. We adopt the G µ scheme with electro-weak parameters taken from the PDG [120], namely The fiducial volume is defined by applying the following set of selection cuts on the lepton pair [97] 8 : where p ± t are the transverse momenta of the leptons, η ± are their pseudo-rapidities in the hadronic centre-of-mass frame, and M is the invariant mass of the di-lepton system. We also define an 'inclusive' setup by dropping in eq. (5.5) the cuts on p ± t and η ± .
Factorisation and renormalisation scales are chosen as µ R = κ R M t , µ F = κ F M t , with M t = M 2 + p t 2 , and p t the di-lepton-system transverse momentum, while the resummation scale is set to Q = κ Q M . For the resummed results, the definition of M t is actually approximated by M , which is appropriate up to quadratic power corrections. We assess the impact of missing higherorder contributions by performing a variation of µ R and µ F by a factor of 2 around their respective central values whilst keeping 1/2 ≤ µ R /µ F ≤ 2. In addition, for central µ R and µ F we vary the resummation scale Q by a factor of 2 in either direction. The final uncertainty for resummed results is built as the envelope of the resulting 9-scale variation, while in the case of matched results, as anticipated above, the envelope also includes variations of the v 0 parameter in eq. (5.2). In Fig. 1 we show a comparison of pure resummed results for the di-lepton transverse-momentum p t distribution in the inclusive setup at NNLL (pink), NNLL without running-coupling effects (orange), NNLL with running-coupling effects (green), and N 3 LL (red). The plot on the left panel displays variations around central scales κ R = κ F = 1, κ Q = 1/2, while the right panel features central scales κ R = κ F = κ Q = 1/2. Both plots clearly show the benefits of the inclusion of 'primed' effects on the NNLL predictions at the level of central value and theoretical-uncertainty bands, especially in terms of shapes. We note that NNLL predictions, both with and without running-coupling effects, are significantly closer to the full N 3 LL result than the NNLL one is, although the pattern of comparison somewhat depends on the chosen central-scale setup, with the running-coupling option closer to full N 3 LL on the left, and the opposite on the right panel. The uncertainty band of the NNLL predictions is also significantly reduced below 10 GeV with    respect to the NNLL one. The band relevant to the running-coupling option is smaller than the the non-running one, which is generally expected since the former encodes correct higher-order running-coupling information, absent in the latter. We note that across the entire p t range the former band is also very similar to the N 3 LL one, and moreover, in all cases does it contain the central N 3 LL prediction, yielding a reliable estimate of the impact of missing higher-order terms. The difference between the two NNLL results may become non negligible at very small p t for certain scale setups (especially so when the central κ Q is different from the central κ R , κ F values), which is also qualitatively expected as due to the approaching to a strong-coupling regime; in all cases the discrepancy is covered by the uncertainty band of the NNLL with running-coupling option, which faithfully assesses the ambiguity related to the inclusion of beyond-accuracy runningcoupling effects. In the following we choose the running-coupling option of 'primed' results as the default for our phenomenological study.
In Fig. 2 we assess the effect of the recoil prescription detailed in Sec. 3.4 on fiducial p t predictions at N 3 LL accuracy (where not explicitly stated, we employ the running-coupling option), both  Figure 3. Left: resummed predictions at N 3 LL (red) and N 3 LL (blue) for p t in the fiducial ATLAS setup. Right: matched prediction at N 3 LL +NNLO (red) and N 3 LL +NNLO (blue). In the right plot, the x axis is linear up to 30 GeV and logarithmic above.
without (left panel) and with (right panel) additive matching (5.1) to the fixed NNLO differential result. The uncertainty band stems from variations around central scales κ R = κ F = 1, κ Q = 1/2, while the matched result includes variation of v 0 as well. The inclusion of recoil (blue, as opposed to purple not featuring recoil effects) gives rise to an expected linear power correction in the pure resummed case, as can be specifically checked in the lower inset of the left panel. After matching to fixed order, recoil induces a marginal ∓1% distortion of the spectrum below 20 GeV, which is the leftover effect after the O(α 3 s ) cancellation taking place between resummation and its expansion in (5.1). The uncertainty bands are also very similar across the whole phase space. Figure 3 displays a comparison, at the fiducial level and including recoil effects, between resummed results (left panel) at N 3 LL (red) and at N 3 LL (blue) accuracy, and between matched results (right panel) at N 3 LL +NNLO (red) and at N 3 LL +NNLO (blue) accuracy. All variations are relevant to central-scale values κ R = κ F = 1, κ Q = 1/2. The inclusion of 'primed' effects on the pure resummed prediction induces a distortion in the spectrum which is less than 2% above 5 GeV, and that can be as large as a few percent below, which is qualitatively consistent with (and quantitatively less pronounced than) what is shown for the NNLL versus NNLL comparison in the left panel of Fig. 1, featuring the same central-scale setup. The uncertainty band undergoes a significant reduction below 10 GeV in passing from N 3 LL to N 3 LL accuracy, by up to a factor of 2 towards p t → 0. The matched results shown on the right panel largely inherit the features just described in the phase-space region dominated by resummation effects, whereas for p t above 50 GeV the prediction is dominated by the fixed-order component, which is common to both. Overall, the N 3 LL +NNLO residual uncertainty band is at the level of 2 -3% below 30 GeV (barring the first bin), and around 5% above 30 GeV. Fig. 4 shows a comparison of the default additive-matching prescription defined in eq. (5.1) (blue) with the multiplicative matching defined in eq. (5.3) (orange) at the N 3 LL +NNLO level, where both predictions include transverse-recoil effects. For reference, the central-scale setup is κ R = κ F = 1, κ Q = 1/2, and the additive prediction is the same as in the right panel of Fig. 3. The theoretical systematics related to the choice of matching family results fairly negligible at this order, with the two predictions being essentially indistinguishable both at central scales, and with respect to uncertainty bands. As the envelope of the two different schemes essentially coincides with the single uncertainty bands, we refrain from adopting it as an estimate of matching systematics,    Figure 5. Comparison of matched predictions at N 3 LL +NNLO (red) and N 3 LL +NNLO (blue) with ATLAS data [97] for p t (left panel) and φ * η (right panel). The fixed-order component is turned off below φ * η = 3.4 · 10 −2 in the right panel, see main text for details. In the left plot, the x axis is linear up to 30 GeV and logarithmic above. and rather insist on the variation of parameter v 0 in a sensible range, such as [2/3, 3/2] around the central v 0 value, as better suited to this aim. This variation is responsible for the slight widening of the band between 30 GeV and 100 GeV, which we believe to reflect a genuine matching uncertainty in this region.
In Fig. 5 we finally compare matched predictions in the fiducial setup to ATLAS data [97], both for p t (left panel) and for φ * η (right panel). The left panel includes the same theoretical predictions shown in the right panel of Fig. 3 (keeping the same colour code), which are here normalised to their cross section in order to match the convention of the shown data. The matched N 3 LL +NNLO predictions for p t show a remarkable agreement with experimental data, with a theoretical-uncertainty band down to the 2 -5% level, essentially overlapping with data in all bins form 0 to 200 GeV (barring one low-p t bin, where the cancellation between the fixed-order and the expanded components is particularly delicate, and few middle-p t bins where the agreement is only marginal). The inclusion of 'primed' effects tends to align the shape of the theoretical prediction to data, so that the former never departs more than 1 -2% from the latter below 200 GeV, as opposed to the more visible relative distortion of the N 3 LL +NNLO below 5 GeV and above 50 GeV. The φ * η results on the right panel follow by and large the same pattern just seen for p t , with 'primed' effects being relevant to improve the data-theory agreement over the entire range, expecially at very small φ * η , and theoretical uncertainties at or below the ±3% level. We incidentally note that, due to the extremely soft and collinear regime probed by φ * η data, the fixed-order component features some fluctuations at small φ * η ; consequently, we have opted to turn it off in the first bins (up to φ * η = 3.4 · 10 −2 ), which implies that the matching formula in that region corresponds to the sole resummation output, multiplied by Z(v). On the one hand this shows that resummation alone is capable of predicting data remarkably well both in shape and in normalisation at very small φ * η ; on the other hand it highlights the necessity of dedicated high-statistics fixed-order runs in order to reliably extract information on fiducial cross sections at N 3 LO by means of slicing techniques, especially in presence of symmetric lepton p ± t cuts.

Higgs results
For Higgs phenomenology we consider hadro-production at the 13 TeV LHC in an inclusive setup, with an un-decayed Higgs boson and no cuts, as well as in a fiducial setup, where we focus on the H → γγ decay channel. We employ an effective-field-theoretical (HEFT) description of the gluon-fusion process where the the top quark running in the loops is integrated out, giving rise to an effective ggH coupling. As seen above, the hard-function coefficients H  where p γi 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 phase-space cuts we do not include the photon-isolation requirement of [121], since this would introduce additional non-global logarithmic corrections in the problem, spoiling the formal accuracy of the resummation. However, we point out that the photon-isolation is quite mild in this particular setup, hence it could faithfully be included at fixed order. The photon decay is predicted in the narrow-width approximation applying a branching ratio of 2.35 × 10 −3 .
For fiducial predictions we employ parton densities from the PDF4LHC15_nnlo_mc set [122][123][124][125][126][127]. Central renormalisation, factorisation, and resummation scales are set as µ R = κ R M H , µ F = κ F M H , Q = κ Q M H , respectively. Theoretical-uncertainty bands are obtained as explained in section 5.1 for the Drell-Yan case. In the inclusive setup, used solely to show the impact of running-coupling effects on 'primed' results, we employ the NNLO NNPDF3.1 PDF set [118] with α s (M Z ) = 0.118.
In Fig. 6 we consider inclusive Higgs production, and show pure resummed predictions for the Higgs transverse momentum p H t at NNLL, NNLL , and N 3 LL with central-scale choices κ R = κ F = 1, κ Q = 1/2 (left panel), and κ R = κ F = κ Q = 1/2 (right panel). This figure, which is the exact analogue of Fig. 1 discussed above, aims at assessing the effect of including or not running-coupling effects in 'primed' results relevant for Higgs production. The benefit of including 'primed' predictions proves significant in this case as well, but with a different pattern with respect to Drell-Yan production. The shape distortion in passing from NNLL to NNLL has a slightly more limited range, mainly extending up to 5 GeV in p H t ; however, the normalisation of the theoretical   Figure 7. Left: resummed predictions at N 3 LL (red) and N 3 LL (blue) for p γγ t in the fiducial ATLAS setup. Right: matched prediction at N 3 LL +NNLO (red) and N 3 LL +NNLO (blue). In the right plot, the x axis is linear up to p γγ t = 50 GeV and logarithmic above.
curves is significantly affected, with 'primed' predictions correctly capturing the large K-factor, at the level of 15% at this perturbative order, which is known to arise in Higgs production. We note the the two different NNLL predictions are fairly similar, and remarkably closer (in shape and normalisation) to the N 3 LL one than the bare NNLL is, both in terms of central value, and of uncertainty-band estimate. The central NNLL prediction without running coupling tends to be slightly closer to the central N 3 LL one, while NNLL with running coupling is slightly more similar to N 3 LL in terms of uncertainty band. In all cases does the central N 3 LL prediction lie well within the NNLL running-coupling band, which we use as our default for the fiducial study. Fig. 7 displays a comparison, relevant to the fiducial di-photon p γγ t spectrum, of N 3 LL curves (blue) agains N 3 LL predictions (red), both without (left panel) and with (right panel) additive matching to NNLO. All predictions include recoil effects, so that this figure represents the Higgsproduction analogue of Fig. 3, but referred to central scales κ R = κ F = κ Q = 1/2. The shape distortion with respect to N 3 LL predictions is more modest in the Higgs case with respect to Drell-Yan production, partly owing to the chosen central-scale setup; moreover, the induced K-factor is fairly close to unity at this order, which is sign of a good perturbative convergence. Overall, N 3 LL predictions feature a significant reduction in theoretical uncertainty in comparison to N 3 LL ones, especially in the low-p γγ t region dominated by resummation. Residual uncertainty is as low as 5 -7% below 10 GeV, and in the matched case it never exceeds 10% below 40 GeV.
Finally, in Figure 8 we show a comparison of theoretical predictions for the fiducial p γγ t spectrum at N 3 LL +NNLO (red) and N 3 LL +NNLO (blue) level, with recoil effects, against ATLAS preliminary data [98]. Theoretical predictions, based on central scales κ R = κ F = κ Q = 1/2, have been rescaled by a factor K rEFT = 1.06584 to account for the exact top-quark mass dependence at LO.

Conclusion
In this article we have presented state-of-the-art differential predictions relevant for colour-singlet hadro-production at the LHC within the RadISH framework, up to N 3 LL +NNLO order. Such a level of accuracy in the resummed component is reached by supplementing the previously available N 3 LL result with the complete set of constant terms of relative order O(α 3 s ) with respect to the Born level. We have documented in detail how the inclusion of such terms is achieved in RadISH, as well as the validation we have performed to confirm the correctness of their numerical implementation. In this article we have focused on neutral Drell-Yan and Higgs production, although we stress that the formalism used here can be straightforwardly applied to the charged Drell-Yan case as well.
We have assessed the behaviour of 'primed' predictions in inclusive Drell-Yan and Higgs production in a comparison of two different NNLL prescriptions (including or not higher-order runningcoupling effects, respectively) with N 3 LL. This has given us confidence on the mutual consistency of the two 'primed' results, and on the reliability of their quoted uncertainty bands, in view of comparing results based on N 3 LL predictions with experimental data. In particular, in all considered cases are the NNLL uncertainty bands capable of encompassing the N 3 LL central prediction, and to correctly estimate the higher-order running-coupling ambiguity underlying the definition of 'primed' accuracy.
The results presented here are fully exclusive with respect to the Born phase-space variables, lending themselves to be flexibly adapted to the fiducial volumes considered in realistic experimental analyses. In order to more accurately simulate the kinematics of the colour-singlet decay products (we considered a lepton pair in the case of Drell-Yan production, and a photon pair in the case of Higgs production), we have consistently encoded in our prediction a prescription for the treatment of the singlet's transverse recoil against soft and collinear QCD initial-state radiation. This includes in our results the full set of linear next-to-leading-power corrections for azimuthally symmetric observables, such as the transverse momentum of the singlet.
The inclusion of transverse-recoil effects, which is performed at the level of differential (as opposed to cumulative) cross sections, and the availability of O(α 3 s ) constant terms in the resummed component, has led to the definition of two differential matching prescriptions, belonging to the additive and multiplicative families, respectively. We have compared the two schemes in Drell-Yan production, and found very good agreement between them, showing that matching systematics are well under control. Variation of matching parameters has anyway been conservatively included in the estimate of the theoretical uncertainties.
Although the ingredients presented above would immediately allow us to quote numbers for the N 3 LO fiducial Drell-Yan and Higgs cross sections, we refrain from doing so in the present article, as in our opinion such a high-precision prediction requires dedicated high-statistics fixedorder runs in order to avoid potential numerical biases. This is especially the case in the context of a slicing technique in presence of symmetric cuts on the transverse momentum of the singlet's decay products. We leave this development for future work.
As a general upshot of the present work, which we have documented both in Drell-Yan and in Higgs production, the inclusion of 'primed' effects is highly beneficial for the stability of the theoretical prediction, leading to a significant reduction in the residual theoretical uncertainty. In the case of our highest-accuracy result, N 3 LL +NNLO, such a reduction can be as large as a factor of 2 in the region dominated by resummation.
For the di-lepton p t spectrum in Drell-Yan production, the N 3 LL +NNLO prediction is shown to improve the comparison to ATLAS data with respect to N 3 LL +NNLO. The theory-data agreement is now at a remarkable level of 1 -2% below 200 GeV, and the residual theory uncertainty is at or below the 2 -5% level in that phase-space region. Same considerations hold for the φ * η observable, with the N 3 LL +NNLO band nicely overlapping with data over the entire range, with leftover uncertainty below ±3%. For Higgs production we find a similar qualitative pattern, with N 3 LL predictions featuring an uncertainty at the level of 5% at very low di-photon p γγ t , and matched N 3 LL +NNLO results well below ±10% accuracy over the entire p γγ t range.
The RadISH code used for the predictions shown in this paper will be made public in due time, and the results are available from the authors upon request. discussions with Carlo Oleari on the axial structure of the Drell-Yan form factor. The work of LR is supported by the Swiss National Science Foundation (SNF) under contract 200020_188464.

A Parton luminosities up to N 3 LL
We report the explicit expression for the parton luminosities employed in the main text, up to N 3 LL accuracy. By defining the coupling factors , which correspond to α s (k t1 ) written in terms of α s (µ R ) at 0, 1, 2, 3 loops, and the standard luminosities can be compactly written as where H (0) (µ R ) = 1, x 1,2 = e ±Y M/ √ s, and Y is the Born-level rapidity of the colour singlet in the centre-of-mass frame of the collision, with energy √ s. As for the luminosities relevant for 'primed' predictions, they assume a different functional form for the running or non-running options, described in the main text. In the running case, we just set L N k LL = L N k+1 LL ; in the non-running case, L N k LL =L N k+1 LL , where theL N k+1 LL corresponds to a luminosity L N k+1 LL with the replacementᾱ s(j) →ᾱ s(j−1) for j > 0.

B V-A structure of the form factor in the neutral Drell-Yan process
In this subsection we discuss the subtleties arising in the extraction of the hard-virtual corrections to the form factor for the neutral-current Drell-Yan process that we will denote with the shortcut Z/γ * , although in the following the dependence upon the final-state leptonic momenta is explicitly taken into account. Before discussing how the hard-virtual coefficients in eq. (3.10) are extracted from the loop corrections to the form factor, it is useful to recall the structure of the Drell-Yan tree-level amplitude expressed in terms of spinor currents. The fermion-antifermion-photon, and fermion-antifermion-Z vertices are defined as respectively, where Q f is the charge of the fermion in units of the positron charge |e|, and θ W being the weak mixing angle and T 3 f = ±1/2 the weak isospin quantum number of the fermion type f . The fermionic currents relevant for the process in eq. (B.1) are where, for the initial-state quark current (f = q), the Dirac spinors readF =v(p 2 ) and F = u(p 1 ), whereas, for the leptonic current (f = ),F =ū(p 3 ) and F = v(p 4 ).
The tree-level amplitude can be written as where q = p 1 + p 2 , and the symbol '·' represents the Lorentz dot product (among currents, in this case). The superscript on the currents indicates the loop order at which they are computed. By denoting the products of currents of fermion line f with the tree-level squared amplitude reads (up to global factors not relevant for the present discussion) where 9  The H (n) (µ R ) coefficients in eq. (3.10) are of pure hard-virtual origin, and they are obtained from the finite parts of the MS-renormalised loop corrections ('r,f' in the equations below) to the tree-level amplitude. For Z/γ * production, if one solely considers loop diagrams where the external quark line is directly connected to the Z/γ * vertex, the full tree-level squared amplitude (B.7) can be factored in front of the hard function, namely up to three loops one has  The radiative corrections to these three current correlators are the same, and they can be extracted from the quark form factor, despite the latter quantity is defined as the coupling of a virtual photon to a massless quark-antiquark pair, and hence, by definition, it contains the radiative corrections to J µ V (q), i.e. to the vector form factor. The physical reason why the radiative corrections to the massless vector and axial-vector form factor are the same ultimately stems from chirality conservation in QCD. Algorithmically, this a consequence of the properties of the γ 5 matrix: the Dirac trace in [J J V V (q)], therefore the A 2 q term receives exactly the same radiative corrections of, e.g., the term Q 2 q [1]. The computation of the radiative corrections to the terms proportional to V q A q and Q q A q is more delicate, as they involve a fermionic trace with just one γ 5 matrix, hence they require a consistent treatment of the axial-vector current in d dimensions, which requires care from two loops [5,128]. Eventually, as happens at tree level, the only nonvanishing contraction of [J J In the light of the above discussion, for Z/γ * production, the tree-level squared amplitude (B.7) can be factored out, as in eq. (B.9), and the H (n) coefficient can be obtained from the quark form factor at n loops.
Starting from NNLO, the vertex corrections to Z/γ * production contain graphs where the external off-shell gauge boson does not couple directly to the external quark line with flavour f = q, but rather to an internal closed fermion loop. In Fig. 9 we show representative diagrams at two and three loops. For such terms, customarily referred to as 'singlet' contributions, the factorisation in (B.9) is violated. Up to N 3 LO, Z/γ * production features two singlet contributions, namely (a) non-vanishing finite corrections entering the axial-vector current but not the vector one, starting from two loops; (b) corrections to the vector and axial-vector currents not factorising the tree-level form factor, at three loops.
As for contribution (a), at two loops the singlet correction to the vector current vanishes identically for each quark running in the fermion loop, by means of Furry's theorem. Since the axial-vector coupling is proportional to T 3 f , within a given generation, the singlet correction to the latter also vanishes, provided the two quarks are degenerate in mass. Such a cancellation does not take place exactly in the third family. The leftover contribution, finite owing to the fact that the Standard Model is anomaly free, and vanishing in the m top → m bottom limit, has been computed in ref. [129]. For a given external quark line q, the radiative correction is proportional to A q q A q J(m q , M Z ) (eq. (7), Ref. [129]). As this is a correction to H (2) that is proportional to [J J (00) AA (q)] but whose coupling is not proportional to A 2 q , it does not factorise |M (0) | 2 . To our knowledge, these contributions are typically not included in resummed calculations. In our implementation, we have not included these axial corrections to the vertex.
We also stress that at O(α 2 s ) there are other terms of this type arising from the 'real-virtual' interference of the + − +1 jet matrix elements. Similarly to the two-loop singlet axial corrections to the vertex, these terms are UV-and IR-finite, and vanish for each mass-degenerate quark family running in the fermion loop. Such 'real-virtual' corrections are also absent from our prediction.
The singlet contribution (b) to the quark form factor at three-loops has been computed in Refs. [77,130,131]. In our notation, it contributes to the third-order expansion of J µ V (q) as α s 2π (B.12) The J (3s) V (q) current couples to the external quark line in the loop amplitude through q Q q (or q V q ) for γ * (or Z) exchange, where index q labels all possible quark flavours running in the closed fermion loop. The singlet contribution enters the hard-virtual coefficient through the interference term 2 (M (0) * M (3s) ) (we drop the subscript r,f as J We note that, in case of γ * production alone, one could collect a term q Q q /Q q in the first line of eq. (B.13), thereby expressing the full result in a factorised form, as done in refs. [77,130,131], which is not possible when considering both γ * and Z channels. By taking the interference with the tree-level amplitude (B.5), one gets where |A (03s) In our implementation we included, at O(α 3 s ) the singlet contribution as in eq. (B.14), with the sum over q running over the 5 massless flavours. We have not included in our calculation the singlet contribution to the massless axial-vector quark form factor at three loops, that has been computed very recently in [132]. We leave this update for a future development, as this contribution is expected to be numerically negligible for the results presented in this article.