Power Corrections for N-Jettiness Subtractions at ${\cal O}(\alpha_s)$

We continue the study of power corrections for $N$-jettiness subtractions by analytically computing the complete next-to-leading power corrections at $\cal{O}(\alpha_s)$ for color-singlet production. This includes all nonlogarithmic terms and all partonic channels for Drell-Yan and gluon-fusion Higgs production. These terms are important to further improve the numerical performance of the subtractions, and to better understand the structure of power corrections beyond their leading logarithms, in particular their universality. We emphasize the importance of computing the power corrections differential in both the invariant mass, $Q$, and rapidity, $Y$, of the color-singlet system, which is necessary to account for the rapidity dependence in the subtractions. This also clarifies apparent disagreements in the literature. Performing a detailed numerical study, we find excellent agreement of our analytic results with a previous numerical extraction.


Introduction
The precision study of the Standard Model at the LHC, as well as increasingly sophisticated searches for physics beyond the Standard Model, require precision predictions for processes in a complicated hadron collider environment. When calculating higher-order QCD corrections, the presence of infrared divergences require techniques to isolate and cancel all divergences. Completely analytic calculations are only possible for some of the simplest cases, e.g. [1][2][3], while for more complicated processes, in particular those involving jets in the final state, numerical techniques are typically required. At next-to-leading order (NLO) the FKS [4,5] and CS [6][7][8] subtraction schemes provide generic subtractions for arbitrary processes, and have been used to great success. At next-tonext-to-leading order (NNLO), due to the more complicated structure of infrared singularities, the development of general subtraction schemes has proven more difficult. While subtraction schemes have been demonstrated both for color-singlet production [9][10][11], as well as for several processes involving jets in the final state [12][13][14][15][16][17][18], significant work is still required before efficient NNLO subtractions can be achieved for arbitrary colored final states.
An important feature of N -jettiness subtractions is that power corrections in the resolution variable can be calculated in an expansion about the soft and collinear limits, allowing the numerical performance of the subtractions to be systematically improved. Recently there has been significant interest in understanding subleading power corrections to collider cross sections [52][53][54][55][56][57][58][59][60][61][62][63][64][65][66][67]. Advances in the understanding of subleading power limits using SCET [61,62] have allowed the leading logarithmic (LL) next-to-leading power (NLP) corrections to be computed at NLO and NNLO [58,60], with independent calculations of the same terms done by a second group in refs. [59,66]. The leading logarithms have also been resummed to all orders for pure glue QCD for 2-jettiness in H → gg [68].
The inclusion of the leading logarithms was found to improve the numerical accuracy, and thereby the computational efficiency, of N -jettiness subtractions for color-singlet production by up to an order of magnitude [58,60]. The analytic calculation of the next-to-leading logarithmic (NLL) power corrections is important for several reasons. Theoretically, it greatly furthers our understanding of the power corrections, since the perturbative structure becomes significantly more nontrivial at NLL as compared to at LL. From a practical perspective, they provide further substantial improvements in the numerical performance of the subtractions. In particular, they make the subtractions much more robust in cases where there are accidental cancellations between different channels or where the NLL terms are numerically enhanced relative to the LL terms.
In this paper, we compute the full NLP corrections at O(α s ) for both Drell-Yan and gluonfusion Higgs production in all partonic channels, including the nonlogarithmic terms (which are the NLL terms at O(α s )). One focus of this paper is to discuss in detail the subleadingpower calculation, in particular the treatment of Born measurements at subleading power, which in our case are the invariant mass Q and rapidity Y of the color-singlet system, and to derive a master formula for the NLP corrections to 0-jettiness at O(α s ). Our analysis lays the ground for extending the calculation of the NLP corrections to higher powers, higher orders in α s , and to more complicated processes.
We also perform a detailed numerical study comparing our analytic results with the full nonsingular result extracted numerically from MCFM8 [23,[69][70][71], which allows us both to verify our analytic calculation and to probe the typical size of the higher-power corrections in various different partonic channels. We find that the NLL power corrections can exhibit a much more pronounced Y dependence from PDF effects than is the case at LL, which demonstrates the importance of calculating the power corrections fully differential in the Born phase space.
Our discussion of the Born measurements also allows us to clarify an apparent disagreement in the recent literature regarding the LL power corrections. As we discuss in more detail in sec. 6, since the calculations in refs. [59,66] are not differential in the color-singlet rapidity Y , their results can only be used fully integrated over Y . 1 In contrast, the results computed here, which agree with those previously derived by a subset of the present authors in refs. [58,60], are differential in Y . When integrating over all Y , integration by parts is used to explicitly show that these LL results are equivalent to those of refs. [59,66]. In sec. 6 we also compare our new differential NLL results with the integrated results of ref. [66].
The outline of this paper is as follows: In sec. 2, we briefly review N -jettiness subtractions for color-singlet production and define our notation. In sec. 3, we discuss in some detail the treatment of Born measurements at subleading power. In sec. 4, we derive a formula to NLP for the soft and collinear power corrections for 0-jettiness. Although our primary focus is on NLO, the general strategy is valid to higher orders as well. In sec. 5, we use our master formula to derive explicit results for the NLP power corrections at NLO for both Drell-Yan and gluon-fusion Higgs production. In sec. 6, we provide a detailed comparison with the literature for those partonic channels where results are available. In sec. 7, we present a detailed numerical study, and compare our analytic results with a previous numerical extraction. We conclude in sec. 8.

N-Jettiness Subtractions, Definitions and Notation
In this section we briefly review N -jettiness subtractions [16,18] in the context of color-singlet production, and discuss the structure of the power corrections to the subtraction scheme. This also allows us to define the notation that will be used in the rest of this paper. For a detailed discussion, we refer the reader to ref. [18].
To compute a cross section for color-singlet production σ(X), where X denotes some set of cuts on the Born phase space, we write the cross section as an integral over the differential cross section in the resolution variable T 0 For a general measure, the 0-jettiness T 0 can be defined as [41,72] where the sum runs over all hadronic momenta k i in the final state. Here, q a,b are projected Born momenta (referred to as label momenta in SCET), which are given in terms of the total leptonic invariant mass Q and rapidity Y as  [19,40] as beam thrust, are the leptonic and hadronic definitions given by leptonic: It has been shown [58] that the power corrections for the hadronic definition are poorly behaved, and grow exponentially with rapidity, while the e ±Y factor in the measure for the leptonic definition exactly avoids this effect. For later convenience, we write the dimensionful and dimensionless 0-jettiness resolution variables in terms of a Y -dependent parameter ρ(Y ) as with leptonic: In the following, we will mostly drop the subscript 0 on T 0 , since there should be no cause for confusion that our results are for 0-jettiness. For generic results that apply to both the leptonic and hadronic definitions we also drop the superscript and simply use T and τ , keeping a generic parameter ρ when necessary.
To implement the N -jettiness subtractions, we now add and subtract a subtraction term to the cross section (suppressing the dependence on the Born measurements X for simplicity) Since T 0 is a zero jet resolution variable, for τ = T 0 /Q → 0 we can expand the differential cross section dσ/dτ and its cumulative σ(τ cut ) about the soft and collinear limits from τ → 0 and τ cut → 0 as Here dσ (0) /dτ and σ (0) (τ cut ) contain all leading-power terms, These terms must be included in the subtraction term to obtain a finite result, namely The further terms in the series expansion in eq. (2.10) are suppressed by powers of τ While these terms with k ≥ 1 do not need to be included in the subtraction term, the size of the neglected term, ∆σ(τ cut ) is determined by the leading-power corrections that are left out of σ sub . Therefore, including additional power corrections in σ sub can significantly improve the performance of the subtraction. Indeed, general scaling arguments imply that up to an order of magnitude in performance can be gained for each subleading power logarithm that is included in the subtractions [18]. For the leading logarithms, this was explicitly confirmed for most partonic channels in the numerical studies in refs. [58,60]. Here, we extend the calculation to the NLL terms at O(α s ), which yields the nonlogarithmic terms, hence giving the complete NLP result. The remaining NLO power corrections then scale at worst as α s τ 2 cut log τ cut , and will be very small, as we will see in our numerical studies.

Born Measurements at Subleading Power
We begin by discussing in some detail the treatment of the Born measurements, Q 2 and Y , which plays an important role at subleading powers. We will use the soft and collinear expansions from SCET, which provide a convenient language when discussing the power expansion of QCD amplitudes at fixed order. We will not need to employ any of the field theory technology from SCET for our analysis here.

General Setup and Notation
Consider the production of a color-singlet final state L of fixed invariant mass Q and rapidity Y , together with an arbitrary measurement T that only acts on hadronic radiation and gives T = 0 at Born level. Since the observable T resolves soft and collinear emissions it will induce large logarithms ln(T /Q). Our goal is to expand the cross section in T (or τ = T /Q) in order to systematically understand its logarithmic structure. Consider proton-proton scattering with the underlying partonic process where L is the leptonic (color-singlet) final state and X denotes additional QCD radiation. Its cross section reads Here, the incoming momenta are given by k = i k i is the total outgoing hadronic momentum, and q is the total leptonic momentum. Since our measurements are not sensitive to the details of the leptonic final state, we have absorbed the leptonic phase space integral into the matrix element, The matrix element M contains the renormalization scale µ 2 , which as always is associated with the renormalized coupling α s (µ), and may also contain virtual corrections. For now the measurement functionT ({k i }) is kept arbitrary.
We can now solve the Q 2 and Y measurements to fix the incoming momenta as Taking the Jacobian factors from solving the δ functions into account, eq. (3.2) becomes dσ where we defined to stress that the squared matrix element only depends on the Born measurements Q and Y , which fix the incoming momenta through eqs. (3.3) and (3.5), and the emission momenta k i . Note that we have left implicit in our notation in eq. (3.6) the dependence of ζ a,b on k through eq. (3.5). They are restricted to ζ a,b ∈ [0, 1], which is implicit in the support of the proton PDFs.

Power Expansion in Soft and Collinear Limits
Instead of solving the T measurement function to express (some of the) k i in terms of T , we find that a convenient strategy to organize the expansion in T is to multipole expand the final state momenta. At this stage we need only assume that T is a SCET I observable, which is true for many definitions of N -jettiness. For such observables, it is known from SCET I that we can organize the cross section in terms of a power counting parameter λ ∼ √ τ . All momenta k i can then be categorized as either collinear or soft modes (since we work in SCET I these are often called ultrasoft, although we will not make this distinction), whose momenta scale as where we decomposed each momentum into lightcone coordinates Here n andn are lightlike vectors satisfying n ·n = 2. The components of the momenta that scale like λ 2 are referred to as residual momenta. The soft momenta are homogeneous, and have purely residual scaling. Overlap between the soft and collinear modes occurring in integrals over final state momenta is removed by the zero-bin subtraction procedure [73]. The benefit of this decomposition is that it allows one to expand eq. (3.6) in λ, agnostic of the actual measurement T . The LP result is then simply obtained by expanding the cross section through λ 0 , the NLP result by expanding through λ 2 , etc. Note that when performing this expansion, all other factors, such as Q, E cm ∼ λ 0 .
While the expansion of the matrix element is of course process dependent, we can give general expressions for the incoming momentum fractions eq. (3.5), independent of the process and observable T . If k is a soft momentum, then the expansion required at NLP is given by where we factored out the Born momentum fractions In the n-collinear limit, we obtain For clarity, we have grouped terms of the same power counting in round brackets. Similarly, one can obtain then-collinear limit, or any combination as might appear when combining multiple emissions.

Master Formula for Power Corrections to Next-to-Leading Power
In this section we derive a master formula for the NLP corrections. This formula applies to any SCET I observable in color-singlet production. In sec. 5, we will apply it to derive explicit results for Drell-Yan and gluon-fusion Higgs production.

General Setup for Color-Singlet SCET I Observables
For reference, we start with the LO cross section for the production of a color-singlet final state L of invariant mass Q 2 and rapidity Y , together with an (up to now arbitrary) measurement T acting only on hadronic radiation, where x a,b = Q Ecm e ±Y and A LO is the squared matrix element in the Born kinematics, see eq. (3.7). For future reference, we also define the LO partonic cross section,σ LO (Q, Y ), by Next, consider an additional real emission to the Born process. Eq. (3.6) yields where we remind the reader that the incoming momenta p a,b are given by eq. (3.5), From these solutions, we see the interesting feature that at subleading power, regardless of the type of final-state emission, the momenta entering both PDFs are modified. Since we do not measure the azimuthal angle of k, it can be integrated over, So far, this expression is still exact. In the next step, we wish to expand the NLO cross section in λ ∼ T /Q. When T is a SCET I observable, we can use the EFT knowledge from SCET I to expand the momentum k in both collinear and soft limits, as discussed in sec. 3. At NLP, we need to expand eq. (4.6) consistently through O(λ 2 ). The O(λ 2 ) power corrections arise from the following sources: • The incoming momenta ζ a,b : While collinear and soft limits yield quite different power expansions, both give a well-defined expansion in λ. We thus simply define the expansion where z a,b ∼ λ 0 and ∆ (2) a,b ∼ λ 2 . We have pulled out the Born momentum fractions x a,b and written 1/z a,b as a fraction for later convenience. Explicit expressions can be obtained from eqs. (3.10) and (3.12), and will be given below.
• PDFs: Since the momenta ζ a,b enter the PDFs, these also have to be power expanded, • Flux factor: Similar to the PDFs, we have to expand the flux factor (4.9) • Matrix element: The expansion of the matrix element depends both on the process and the considered limit. Here, we define the LP and NLP expansions by (4.10) In the soft limit A (0) ∼ λ −4 and A (2) ∼ λ −2 , while in the collinear limit A (0) ∼ λ −2 and A (2) ∼ λ 0 . In both cases we have the scaling dk + dk − A (2j) ∼ λ 2j , which is why the soft and collinear corrections enter at the same order.
For example, for a 2 → 2 process, the matrix element can be written in terms of the Mandelstam variables Since these terms now have a definite power counting, one can simply insert eq. (4.11) into A(Q, Y ; {k}) and expand to the required order in λ.
• Measurement: Depending on the observable, the measurement functionT may also receive power corrections. Since 0-jettiness is defined in terms of n,n, Q and Y , none of which receive power corrections in our approach, we do not have such corrections, and will therefore not write them explicitly in the following formulae. More generally, these could be obtained from expanding δ T −T ({k}) .
Inserting all these expansions into eq. (4.6) and expanding consistently to O(λ 2 ), we obtain the LP result (4.12) and the NLP master formula a,b (k) are defined by eq. (4.7). Note that the LP limits of the matrix elements are universal, and hence eq. (4.12) holds independently of the process, i.e. it only depends on the observable T . Although the focus of this paper is on the power corrections, in appendix A we provide a brief derivation of the leading-power terms. Likewise, the A (0) term together with the square bracketed factor on the second line of eq. (4.13) is universal, such that all the process dependence arises from the last A (2) term. We will discuss this in more detail in sec. 4.4.
In the following, we will evaluate eq. (4.13) in both the soft and collinear limit for 0jettiness, eq. (2.7), whose measurement function for one emission is given by (4.14) The value of ρ depends on the specific definition of T , as given in eq. (2.8).

Collinear Master Formula for 0-Jettiness
The expansion of the incoming momenta ζ a,b for an n-collinear emission k ∼ (λ 2 , 1, λ) has been given in eq. (3.12), so the explicit expressions for the expansion eq. (4.7) are Since an n-collinear emission satisfies k − k + , the 0-jettiness measurement eq. (4.14) simplifies to Note that the integration in eq. (4.13) also includes the region k − → 0, where the assumption k − k + is invalid. Indeed, this region corresponds to the soft expansion. It is guaranteed by the zero-bin subtraction procedure that this overlap regime between the soft and collinear limits is not double counted [73]. An important benefit of 0-jettiness and our setup here is that the zero-bin contribution that removes the overlap is scaleless and vanishes in pure dimensional regularization, such that we do not need to consider it further.
Eq. 4.17 fixes the k + integral in eq. (4.13). It is also useful to write the remaining k − integration in terms of z a using eq. (4.16), giving Here the lower bound on the integration follows from the physical support of the PDF f a (x a /z a ). Plugging back into eq. (4.13), we obtain the n-collinear master formula where k is given by It only remains to plug in the expansions of the matrix element A (0) and A (2) and to expand in . Note that even in the n-collinear case considered here, both the n andn-collinear incoming momenta receive power corrections, see eq. (4.4), leading to derivatives of both PDFs in eq. (4.19). The analogous results for then-collinear limit are obtained in the same manner, giving where k is given by

Soft Master Formula for 0-Jettiness
The expansion of the incoming momenta ζ a,b for a soft emission k ∼ (λ 2 , λ 2 , λ 2 ) has been given in eq. (3.10), so the explicit expressions for the expansion eq. (4.7) are Plugging back into eq. (4.13), we get Here, the measurement is given by eq. (4.14), We can further simplify eq. (4.25) by utilizing the fact that A (0) and A (2) have a well defined dependence on k + or k − because of power counting and mass dimension, Here the A's are process-dependent expressions that depend on the Born measurements Q and Y , but are independent of both k + and k − . This implies that the k ± integrals in eq. (4.25) have the generic structure We then find the soft NLP master formula

Universality of Power Corrections for 0-Jettiness
Having derived our master formulas, in this section we comment on the universality of the power corrections. In both the collinear and soft limits the power corrections arising from the derivatives of the PDFs and from the expansion of the flux factor are proportional to the LP matrix element A (0) (Q, Y ), see eqs. (4.19) and (4.29). Since the factorization properties of A (0) (Q, Y ) are universal, most of the NLP corrections are universal as well, in the sense that they essentially reduce to the LO cross section times a universal factor, as we will make explicit below. The only process-dependent piece arises from the NLP expansion A (2) (Q, Y ) of the matrix element. We stress that this limit is defined in our particular choice of Born measurements Q 2 and Y . Using different observables, e.g. q ± , the NLP corrections to ζ a,b in eq. (3.5) would change, inducing also a change of the NLP matrix element.

Universality of Collinear Limit
We begin by considering the n-collinear limit of a real emission amplitude in detail. We consider the Born process where κ i denotes all quantum numbers, including flavor, of the incoming partons, and L is the leptonic final state of momentum q = q a + q b . The incoming momenta for the hard collision are given by Now consider that parton a arises from an n-collinear splitting of a parton with flavor a , To describe this at leading power, we only need the O(λ 0 ) relations for the momenta of the incoming partons, which can be read off from eqs. (4.15) and (4.16), The n-collinear emission is given by eq. (4.20), It follows that the leptonic momentum q = q a + q b − k = q + O(λ 2 ) is equal to the Born momentum q = q a + q b , and hence the collinear splitting does not affect the leptonic phase space at LP. The LP limit only exists if the splitting κ a → κ a + κ 1 is allowed, in which case it is given by the O(λ −2 ) piece of the squared amplitude, where the 1/k + gives rise to the λ −2 behavior of the amplitude. Here the P aa are the -dimensional splitting functions which are summarized in appendix A.
Recall that in our case, the measurement fixes k + = T /ρ. In the notation of eq. (4.10), we hence have for the LP matrix element These results enable us to explicitly give the universal part of the NLP result in the collinear limit. Inserting into the collinear master formula eq. (4.19) and converting to the MS scheme, we find Here, we factored out the LO partonic cross section eq. (4.2), which is only possible because the collinear splitting leaves the leptonic momentum invariant at LP. We have made explicit the universal piece and nonuniversal components. As already discussed, the full nonuniversal structure arises from the NLP matrix element A (2) (Q, Y, {k}) in the last line. It would also be interesting to understand if there is a universal structure to A (2) (Q, Y ). This has recently been studied in ref. [74] for pure n gluon scattering amplitudes at the level of the Cachazo-He-Yuan scattering equations [75,76], where it was proven that in the subleading power collinear limits the tree-level amplitude factorizes into a convolution of the n − 1 gluon integrand and a universal collinear kernel. It would be interesting to understand this at the level of the amplitude itself, as well as for fermions. Unlike at leading power, we do not expect that there are universal subleading power splitting functions that are simply functions of z, but there may exist splitting functions that involve differential or integral operators, as occurs in the soft limit at subleading power [77,78]. Understanding this will be particularly important for generalizing the calculation of the power corrections to more complicated processes.

Universality of Soft Limit
As for the collinear case, the LP soft limit of the matrix element is universal. Following similar steps as in sec. 4.4.1, one can express the LP soft limit by which only exists for ab = gg, qq and where C = C A , C F is the appropriate Casimir constant. We thus obtain (4.39) As for the collinear case, this emphasizes that the terms arising from the expansion of the PDFs and flux factor are universal, in the sense that they only depend on the universal LP soft limit of the amplitude. The only nonuniversal contributions are |A (2) ± | 2 . However, these terms can in fact be derived from universal formulae [77][78][79] involving differential operators. This has been recently studied in the threshold limit, where one only requires soft contributions [63]. However, when one is away from the threshold limit as considered here, one in general requires collinear contributions, which as discussed above, are not (yet) known to be universal.

Power Corrections at NLO for Color Singlet Production
In this section we give explicit results for the full NLP correction for 0-jettiness at NLO for Higgs and Drell-Yan production in all partonic channels. Since we only consider cases that are s-channel processes at Born level, the LO matrix element only depends on Q and one can factor out the LO partonic cross sectionσ LO (Q). We write the NLP cross section as where as always We will always express the real emission amplitudes in terms of the Mandelstam variables This allows us to straightforwardly obtain the LP and NLP expansion using eq. (4.11). We will give an explicit example of the derivation of the soft and collinear master formulas for the gg → Hg channel, and only summarize the results in the other channels.

Gluon-Fusion Higgs Production
We begin by considering Higgs production in gluon fusion in the m t → ∞ limit. At NLP, there are three different partonic channels, gg → Hg, qq → Hg and qg → Hq, which we consider separately. The calculation for gg → Hg is shown in full detail as an illustration of our master formulae. The LL power corrections were computed in [59,60]. Ref. [60] also computed the qq → Hg NLL power corrections. The NLL power corrections for all partonic channels for gluon fusion Higgs were computed in [66]. We will compare with these results in sec. 6. Throughout this section we consider on-shell Higgs production, for which the partonic cross section is given bŷ The LO matrix element in d = 4 − 2 dimensions is given by [80,81] The spin-and color-averaged squared amplitude for g( n-Collinear Limit Expanding eq. (5.6) using eqs. (4.11) and (4.16), the LP and NLP limits of the matrix element are obtained as Since our scaling variable is λ ∼ T /Q, we clearly see that A (0) ∼ λ −2 and A (2) ∼ λ 0 , as required at LP and NLP. Inserting these expansions into eq. (4.19) and converting to the MS scheme yields To expand this in , we collect all powers of (1 − z a ) and then use the distributional identity We also combine the two separate f a f b pieces, as at this level there is no use to further distinguish the universal and process dependent pieces. This yields (5.10) Comparing to eq. (5.1), we can read off the n-collinear kernels, Soft Limit To expand the matrix element in the soft limit, we use eqs. (4.11) and (4.24) to obtain Note that the first term scales as (k The NLP term in the expansion of the amplitude thus vanishes, and in the notation of eq. (4.27) we have Inserting into eq. (4.29) and converting to the MS scheme yields Since there is no NLP matrix element, one can also obtain this from the universal expression for the soft limit in eq. (4.39). Comparing to eq. (5.1), we can read off the soft kernel, Final Result Adding the n-collinear kernel eq. (5.11), then-collinear kernel which follows from symmetry, and the soft kernel eq. (5.15), all poles in cancel as expected, and we obtain Substituting these results into eq. (5.1) yields the NLP cross section for gg → Hg at NLO.

gq → Hq
The gq → Hq channel has power corrections at both LL and NLL. The spin-and coloraveraged squared amplitude for g(p a ) + q(p b ) → H(q) + q(k) is given by [80] A Soft Limit The LP soft limit vanishes, since a leading-power soft interaction (which is eikonal) cannot change a n-collinear quark into a n-collinear gluon and soft quark. However this does occur at NLP in the soft expansion and yields 18) and the soft kernel is given by n-Collinear Limit Then-collinear limit has both a LP and NLP contribution, given by (5.20) Then-collinear kernel is obtained as n-Collinear Limit The n-collinear limit vanishes at LP. The NLP expansion of the matrix element gives The only nonvanishing kernel is

(5.23)
Final Result Combining the n-collinear,n-collinear, and soft kernels, the 1/ pole vanishes, and we obtain the final results, Substituting these results into eq. (5.1) yields the NLP cross section for gq → Hq at NLO.

qg → Hq
The final results needed in eq. (5.1) for qg → Hq follow from eq. (5.24) by flipping a ↔ b, e Y /ρ ↔ ρ/e Y and f g ↔ f q ,

qq → Hg
The qq → Hg channel first contributes at NLL. It was first given in [60] and then in [66], which agreed, but we reproduce it here for completeness. The squared matrix element, including the average on the initial state spin and colors, is given by [80] A With our choice of Born measurements, the soft limit vanishes both at LP and NLP, leaving only the collinear NLP correction. The LP collinear limit also vanishes, leaving only the NLP n-collinear limit and then-collinear result is obtained by replacing z a ↔ z b . Combining both, we obtain the kernel for eq. (5.1) -20 -

Drell-Yan Production
We now consider the Drell-Yan process pp → Z/γ * → l + l − , and for brevity denote it as pp → V . At NLO we have the partonic channels qq → V g and qg → V q. The LL power corrections for these channels were calculated to NNLO in [58,59]. For Drell-Yan, it is important to be able to include off-shell effects. The LO partonic cross section as a function of the leptonic invariant mass Q is given bŷ Here, v l,q and a l,q are the standard vector and axial couplings of the leptons and quarks to the Z boson, and we have integrated over the l + l − phase space.

qq → V g
We first consider the partonic channel qq → V g. The squared amplitude is given by [82] |M Soft Limit With our setup, the soft limit of the matrix element has no NLP correction, and the soft kernels are given by Collinear Limit The n-collinear expansion of the matrix element yields (at NLP, we only need → 0) The n-collinear kernel is Final Result Adding the n,n and s kernel, we get Substituting these results into eq. (5.1) yields the NLP cross section for qq → V g at NLO.

qg → V q
Next we consider the partonic channel qg → V q. The squared amplitude is given by [82] A qg→V q (Q, (5.37) Soft Limit The LP soft limit vanishes, and the NLP soft expansion is given by The soft kernel is given by n-Collinear Limit The n-collinear limit does not contribute at LP, since the LP interaction can not change then-collinear gluon into an-collinear antiquark. The NLP matrix element is given by 40) and the collinear kernel is n-Collinear Limit Then-collinear limit is IR finite, so we work in d = 4, Then-collinear kernel is given by Final Result Adding the s, n,n kernels, the pole in cancels and we get Substituting these results into eq. (5.1) yields the NLP cross section for qg → V q at NLO.

gq → V q
For completeness, we also give the explicit results for the gq → V q channel, which can easily be obtained from eq. (5.45) by flipping a ↔ b, e Y /ρ ↔ ρ/e Y and f q ↔ f g ,

Comparison with Integrated Results in the Literature
In this section, we compare our NLO results to previous results in the literature. The LL results presented by a subset of the present authors in refs. [58,60] fully agree with the results obtained in this paper.
The results in refs. [59,66] are given only integrated over the color-singlet rapidity Y , and hence take quite a different form at the integrand level. To compare to them, we integrate our results over Y , which allows us to use integration by parts to bring our results into the same integrated form as those in refs. [59,66]. For leptonic T , whose definition involves Y , we find that ref. [66] uses a different definition, and hence we cannot make a meaningful comparison. For hadronic T , whose definition is independent of Y , we find explicit agreement for the LL results after integrating over Y .
At NLL, the results obtained here for the power corrections differential in Y , for both the leptonic and hadronic definitions and all partonic channels, are new. After integrating over Y we find almost complete agreement with the hadronic results of ref. [66], up to a relatively simple term.
Since there are a number of differences in our treatment compared to refs. [59,66], we provide a detailed comparison in this section. In sec. 6.1 we discuss our different treatments of the NLO phase space and of the Born measurements, and show that the rapidity dependence cannot be easily reconstructed from the results in refs. [59,66]. In sec. 6.2 we provide an explicit comparison of the results for the gg → Hg channel integrated over rapidity at LL and NLL, both analytically and numerically.

Treatment of the NLO Phase Space
The derivation in ref. [66] differs from ours here (and that in refs. [58,60]) in that it is not differential in the rapidity Y . To explore the differences arising from this, we give a brief derivation of the NLO phase space following the same steps as ref. [66]. Note that in the following we always work with an on-shell process, in contrast to our more general setup in sec. 4. We also only consider the case k + < k − , since the case k + > k − follows by symmetry.
We start with the expression for the NLO phase space as given in ref. [66], where s = E 2 cm , Q a is defined in eq. (2.6) and x a arises from the T measurement. We can derive a similar expression in our notation, including in addition the rapidity measurement as done in our main derivation. Denoting the incoming momenta at NLO by q a,b , we have from eq. (3.2) As in eq. (6.1), we assume that k + < k − to setT (k) = ρk + , which gives

Following ref. [66], we now change variables via
Up to the rapidity measurement from the final δ function, we find complete agreement with eq. (6.1) if we identify At this step, our treatment differs from the one in ref. [66]. Since we explicitly implement measurement δ functions for both Q and Y , we can uniquely solve for ξ a and ξ b in terms of Q and Y or equivalently x a and x b , This holds for both ρ = 1 and ρ = e Y . This is equivalent to eq. (3.5) (where we used the notation ζ a,b instead of ξ a,b here). The reason this expression looks different is just because in eq. (3.5) we performed this step before fixing k + in terms of T and before changing variables from k − to z a via k − = ξ a E cm (1 − z a ). Following a similar strategy as in sec. 4, one can now replace ξ a,b in eq. (6.4) by the solution eq. (6.6), take the Jacobian from solving the δ functions into account, and then simply expand in T . The main difference to the derivation in sec. 4 is that here, one directly expands the phase space in T , while in sec. 4 we expanded in terms of the generic powercounting parameter λ.
In ref. [66], there is only the Q 2 measurement but no rapidity measurement, i.e. Y is implicitly integrated over. Hence, there is only one constraint for the two variables ξ a , ξ b , whose solution is not unique. They choose to perform the variable transformation from ξ a,b to new variablesx a,b defined by We writex a,b here to distinguish these from the Born variables x a,b = Qe ±Y /E cm that appear in the Born-projected momenta in eq. (2.4). While they satisfyx axb = Q 2 /E 2 cm due to the Q 2 measurement constraint, (1/2) ln(x a /x b ) is not equal to the rapidity Y , which would require the solution in eq. (6.6).
In ref. [66], thex a,b defined by eq. (6.7) also enter in the definition of the 0-jettiness measure in eq. (2.3) in place of x a,b . As a result, the leptonic T definition in ref. [66] is not the same as the usual one with ρ = e Y that we use. Their hadronic definition is the same as ours, as it has no x a,b dependence. Therefore in the following we restrict our comparison to the hadronic definition.
We also note that one cannot easily recover the rapidity dependence from the integrands of the final results in ref. [66]. To see this explicitly, consider inserting the rapidity measurement by comparing eqs. (6.1) and (6.4), which gives On the right-hand side we have carried out the power expansion about T → 0 and the superscript (0) denotes the results for these variables at LP, while ξ (0) a = dξ a /dT T →0 , etc. This accounts for the fact that in general the ξ a,b can depend on T themselves. Equation (6.8) shows that one cannot use the LP expression δ ] to recover the Y dependence from thex a,b dependence of the results in ref. [66], as this does not account for the additional power corrections induced by the Y measurement in the second line of eq. (6.8). This implies that the results in ref. [66] and also those in ref. [59] cannot be used when being differential in rapidity or integrated over bins of rapidity, but only integrated over all Y . This was also confirmed to us by the authors.

Explicit Comparison to Results in the Literature for gg → Hg
Our final results take a quite different form than those in refs. [59,66]. For us, both ξ a and ξ b receive power corrections resulting in derivatives for both PDFs. In contrast, the variable transformation in eq. (6.7) for the case of k + < k − does not yield power corrections for ξ b and hence no derivatives of f b , while the expansion of ξ a yields derivatives of f a (and vice versa for k + > k − ). Due to this different form, one cannot directly compare the integrands of the two results, but one needs to use integration by parts to bring the results into the same form, as we will now show explicitly. In particular, we will show that the results of refs. [58,60], obtained also here, do agree with the results of refs. [59,66] at LL when integrating over all Y .
Integrating our result over Y , and transforming the integration variables to x a,b = Qe ±Y /E cm , we obtain from eq. (5.1) We will show the integration by parts explicitly for the f i f j piece. Let us denote the piece we wish to integrate by parts by D (2,1) , which can be chosen freely. To integrate over Y 1 < Y < Y 2 , we switch the integration variables x a , x b back to Q 2 and Y , use that and integrate by parts with respect to Y . Combining the resulting pieces with those in eq. (6.9), we find The dependence on D (2,1) exactly cancels in this expression. We can choose D (2,1) freely to obtain different forms of the Y -integrated result. The last term in eq. (6.11) is the boundary contribution, which vanishes as Y 1,2 → ±∞, i.e. only if one is fully inclusive in Y . They do in general contribute when placing acceptance cuts on Y .
We now work out explicitly the required integration by parts both at LL and NLL to bring our results into the integrated form as given in refs. [59,66]. For concreteness, we focus on the gg → Hg channel. For the reasons mentioned earlier, we can only compare the results for the hadronic T definition.

Comparison at LL
At LL, our results in eq. (5.16) simplify to These agree with the earlier results obtained by a subset of the current authors in refs. [58,60]. Note that for strict LL accuracy, one can also write the logarithms as ln(T /Q) ± ln(ρ/e Y ) and only keep the ln(T /Q) at LL, while including the ± ln(ρ/e Y ) pieces in the NLL contributions.
(This is the convention used in refs. [58,60] and in sec. 7.) Here, we keep them as part of the LL result, as they are relevant for the comparison with ref. [66]. Up to a trivial change in notation, the LL result given in ref. [66] for hadronic T is As discussed before, thex a,b here are not equal to the Born variables x a,b . Inserting our LL result in eq. (6.12) with ρ = 1 into eq. (6.9), we have where e Y = x a /x b . At the integrand level, the two results clearly have a different form, as was also remarked in refs. [60,66].
To show explicitly that eqs. (6.13) and (6.14) do agree, we integrate by parts to move the f g f g contribution in eq. (6.14) into the f g f g and f g f g terms. Using eq. (6.11), we can achieve this by choosing terms. In both cases, our result in eq. (6.14) and the result of ref. [66] in eq. (6.13) agree. The small difference in the second case arises due to the fact that e ±Y /m H is not exactly the same asx a,b E cm .
Integrating over Y 1 < Y < Y 2 and using that e Y = x a /x b and m H = √ x a x b E cm , eq. (6.14) becomes dσ (2,1) The first two lines exactly reproduce eq. (6.13). The following two lines are the boundary term from integration by parts, which vanishes as Y 1,2 → ±∞. The last line is a NLL effect and can be neglected for the LL comparison. (It is induced by the integration by parts acting on the Y dependence kept inside the argument of the logarithms.) Therefore, the two expressions in eqs. (6.13) and (6.14) agree at LL and at integrated level if and only if one integrates over all rapidity.
To illustrate this numerically, the Y -integrated results are compared in fig. 1. First note that the hadronic LL results in eq. (6.14) do not exactly correspond to those previously given in refs. [58,60]. This is due to the formally NLL terms proportional to ln(ρ/e Y ), discussed below eq. (6.12), which are dropped in the strict LL results in refs. [58,60], but are kept in eq. (6.14). The analogous NLL terms proportional to ln(x a,b E cm ) are also kept in refs. [59,66] and eq. (6.13). Dropping these NLL terms in eqs. (6.13) and (6.14), our and their LL results defined in terms of the same ln(T /m H ) agree exactly, as shown by the solid red and blue dashed curves in fig. 1. 2 The long-dashed orange and dotted blue curves in fig. 1 show the results when using instead ln(T e ±Y /m H ) or ln[T /(x a,b E cm )] to multiply the LL coefficients. The observed difference to the solid red/dashed blue strict LL result has the size of a typical NLL contribution. There is also a very small difference between the long-dashed orange and dotted blue results due to the fact that e ±Y /m H is not exactly the same asx a,b E cm . This difference is exactly accounted for by the last line in eq. (6.16).

Comparison at NLL
We now extend our comparison of the Y -integrated results to NLL, focusing again only on the gg → Hg channel, which contains all possible complications. The full NLL result of ref. [66] can be written as dσ NLP [66] dT cm = To bring our result into this same form, we need to integrate by parts twice, first with respect to Y as shown in eq. (6.11), and then with respect to z a . The details of this calculation are given in appendix B. The final result is shown in eq. (B.5) and is given by the result of ref. [66] in eq. (6.17) plus an extra contribution, dσ (2,1) dT cm = dσ NLP [66] dT cm (6.18) The two results should agree exactly upon integration, and we have not been able to find a source for this discrepancy. As discussed in more detail in sec. 7, the numerical comparison with MCFM provides a strong confirmation of our result. The numerical extraction of the integrated NLL coefficient yields −0.460±0.026, which agrees well with our analytic predicted value of −0.466 (see table 2 below). Dropping the term in the final line of eq. (6.18) would instead predict the value −1.669. 3

Numerical Results
In this section we study our results numerically, including the size of the power corrections and the rapidity dependence. We also compare our analytic results for the O(α s ) NLP power corrections with the full nonsingular spectrum obtained numerically from the LO V +jet and H+jet calculations in MCFM8 [23,[69][70][71]. In refs. [58,60], the NLP corrections were extracted numerically by using a fit of the known form of their logarithmic structure to the nonsingular spectrum from MCFM8. In refs. [58,60], these fits were carried out for the leptonic definition. Here, we have in addition performed the fits also for the hadronic definition. We find excellent agreement between the analytically predicted values and the numerically extracted values for all coefficients, i.e., for the LL and NLL coefficients in all partonic channels for both the leptonic and hadronic definition. This provides a strong and independent cross check for the correctness of the analytic NLL results obtained here. By comparing the complete nonsingular spectrum with our NLP result, we can also assess the importance of power corrections beyond NLP. The NLO power corrections for each partonic channel are extracted from the nonsingular spectrum by using the fit function with τ ≡ T 0 /m Z for Z production and τ ≡ T 0 /m H for Higgs production. Details of the fitting procedure have been described already in refs. [58,60], so we do not repeat them here. A key point is that in order to obtain a precise and unbiased fit result for the to-be extracted a i coefficients, it is crucial to include the higher-power b i and c i terms in eq. (7.1), and to carefully choose the fit range and verify the stability of the fit, as was done in refs. [58,60]. At the level of precision the a i are extracted, this is essential since the full nonsingular cross section includes the complete set of power corrections and if the b i and c i terms were neglected, these higher-power corrections would be absorbed by the a i terms in the fit, rendering their numerically extracted values meaningless. To obtain a precise extraction of the NLL coefficient a 0 , we fix the LL coefficient a 1 in the fit to its analytic result. The relevant coefficients for our NLP comparison at NLO are the LL coefficient a 1 and the NLL coefficient a 0 . For leptonic T they were extracted for Drell-Yan in ref. [58] and for gluon-fusion Higgs in ref. [60] and for the hadronic T we have obtained them here. Depending on the partonic channel, the uncertainties on the fitted coefficients range from 0.08% to 2.3% for leptonic T and from 0.6% to 5.7% for hadronic T . The latter has larger uncertainties because its power corrections are larger, requiring the fit to be restricted to smaller T values where the uncertainties in the nonsingular data are larger.

Drell-Yan Production
We first consider Drell-Yan production, taking pp → Z/γ * at E cm = 13 TeV. We use the MMHT2014 NNLO PDFs [83] with fixed scales µ r = µ f = m Z , and α s (m Z ) = 0.118. We fix Q = m Z , integrate over the vector-boson rapidity, and work in the narrow-width approximation for the Z-boson. The NLP corrections for the leptonic T definition were numerically extracted in ref. [58]. The results for both the leptonic and hadronic definitions for all partonic channels are collected and compared to our analytic predictions in table 1. We find excellent agreement within the fit uncertainties in all cases.  Table 1. Comparison between our analytic predictions and the fitted results for the LL a 1 and NLL a 0 coefficients in Drell-Yan production. These fitted values for a 1 and a 0 with the leptonic definition and the analytic results for a 1 were already given in ref. [58].
In fig. 2 we show the complete NLO nonsingular contributions as black dots, as well as a fit to their form with the solid red curve. Given the agreement in table 1 between our analytic a 0 and the earlier fit result for a 0 , we have fixed a 0 to the analytic result, and redone the fit using eq. (7.1) to obtain this red curve. The red curve from this fit is fully consistent with the earlier fit result from ref. [58]. The dashed orange curve in fig. 2 is the extension of the fit function beyond its fit range. In dotted green and dashed blue we show our analytic predictions. We see that with the inclusion of the NLL power corrections, we obtain an excellent description of the full nonsingular cross section up to nearly T 0 ∼ 1 GeV. This is quite remarkable, and shows that additional higher-order power correction terms are truly suppressed. In fig. 3 we show a plot of the corresponding residual power corrections for the cumulant, ∆σ(τ cut ), on both a linear scale (left) and logarithm scale (right). The solid red curve shows the full power corrections, the solid green curves show the remaining power corrections after including a 1 in the subtractions, and the solid blue curve those after including a 1 and a 0 in the subtractions. We see that with the inclusion of the full NLL power corrections, we achieve more than a factor of 100 reduction in the residual power corrections as compared with the leading-power result at NLO. Both partonic channels have similarly sized power corrections and show a fast convergence of the power expansion. The fact that the blue curve in the logarithmic plot exhibits a steeper slope than the red and green curves is due to its O(τ 2 cut ) scaling corresponding to a next-to-next-to-leading power correction. This provides a nice visualization that our results correctly capture the complete NLP contribution.
The analogous results for the fitted nonsingular spectrum and the residual power corrections ∆σ(τ cm cut ) for the hadronic T definition are shown in figs. 4 and 5. As expected, the power corrections are substantially larger for T cm than for the leptonic definition. To obtain similarly sized power corrections, one has to go to about an order of magnitude smaller values of T cm . Apart from the overall enhancement, the qualitative behavior of the LL and NLL contributions and the different partonic channels is the same. This is expected from our analytic results, which show that the coefficients for both definitions have essentially the same structure and primarily differ in the overall factors of e ±Y leading to the rapidity enhancement for the hadronic definition already observed in refs. [58,60].
In fig. 6 we show the rapidity dependence of the NLP corrections at fixed τ cut = 10 −3 for both leptonic and hadronic T normalized to the LO rapidity spectrum. We can clearly see the exponential enhancement for the hadronic definition at large |Y |. For the qg channel, the asymmetric behavior in rapidity is expected from its analytic result. The result for the gq channel corresponds to taking Y → −Y , such that their sum is symmetric in rapidity. While the leptonic definition does not suffer from the exponential enhancement of the hadronic definition, it still exhibits a substantial increase at large positive Y in the qg channel, as well as a suppression at large negative Y . This is due to the substantially different x-dependence of the quark-gluon luminosity (and its derivative) compared to the qq luminosity in the LO result to which we normalize. Knowing the NLL contribution to the power corrections differential in rapidity enables one to explicitly account for this effect in the subtractions.   Figure 6. The NLO NLP corrections as a function of rapidity at fixed τ cut = 10 −3 for Z production for the qq channel (top row) and the qg channel (bottom row). The LL and NLL coefficients for leptonic T are shown by the green dotted and blue dashed curves and for hadronic T by the dotted and dashed gray curves.

Gluon-Fusion Higgs Production
Next, we consider gluon-fusion Higgs production. We take pp → H at E cm = 13 TeV with an on-shell, stable Higgs boson with m H = 125 GeV, integrated over all Y . We use the MMHT2014 NNLO PDFs [83], with fixed scales µ r = µ f = m H , and α s (m H ) = 0.1126428. The NLP power corrections for this configuration for the leptonic T definition were extracted numerically in ref. [60]. The results for both leptonic and hadronic definitions for all partonic channels are collected and compared to our analytic predictions in  Table 2. Comparison between our analytic predictions and the fitted results for the LL a 1 and NLL a 0 coefficients in Higgs production. These fitted values for a 1 and a 0 with the leptonic definition and the analytic results for a 1 were already given in ref. [60].
In fig. 7 we show as the solid red curve a fit to the full nonsingular result at NLO (black points), which is compared with the LL and NLL predictions in dashed green and dashed blue, respectively. Once again this solid red fit curve is obtained using the form in eq. (7.1) with a 1 and a 0 fixed by the analytic result in table 2, and agrees very well with the corresponding result obtained in ref. [60] where a 0 was a parameter in the fit. In all cases, we find that the NLL result provides a good description of the full nonsingular cross section. This is expected since the NLL results includes all NLP terms in the NLO cross section. We see, however, that particularly for the gq + qg channel, the NLL result for a 0 is required to get a good description, and the LL power correction a 1 alone is not sufficient. Thus the gq + qg channel provides an example where simply looking at the size of the residual nonsingular result after subtracting the a 1 term does not suffice to validate the value of this coefficient.
In fig. 8, we show a plot of the corresponding power corrections for the cumulant, ∆σ(τ cut ), on both a linear scale (left) and logarithm scale (right). Here we more easily see that the LL coefficient is numerically suppressed, while in contrast its NLL coefficient is quite larger. Due to this unusual behavior, the NLL result is required to consistently reduce the power corrections as compared with the leading-power result. In the qq channel there is no a 1 term, and significant improvement is apparent from including a 0 .
The analogous results for the fitted nonsingular spectrum and the residual power correc-  Figure 11. The NLO NLP corrections as a function of rapidity at fixed τ cut = 10 −3 for Higgs production for the gg channel (top row) and the gq channel (bottom row). The LL and NLL coefficients for leptonic T are shown by the green dotted and blue dashed curves and for hadronic T by the dotted and dashed gray curves. tions ∆σ(τ cm cut ) for the hadronic T definition are shown in figs. 9 and 10. The power corrections are noticeably larger, though the effect of the rapidity enhancement is not as pronounced as for Drell-Yan, since here the PDFs suppress the cross section contributions at larger rapidities. For the dominant gg → Hg channel there are also numerical cancellations in the NLL coefficient. More precisely the value for a 0 in table 2 arises as a 0 = 2.356+(−2.822) = −0.466, where the first term corresponds to the rapidity-enhanced version of the leptonic a 0 while the second term is the NLL contribution arising from the additional rapidity dependence in the argument of the leading logarithm discussed below eq. (6.12). As a result of this cancellation, including only a 1 in the subtractions leads to slightly smaller power correction above τ cut > 10 −3 than subtracting both a 1 and a 0 (compare the green and blue solid lines in the top row of fig. 10). If the second NLL contribution were included as part of the LL result, the latter would provide a much poorer approximation and including the remaining NLL contribution would provide a substantial improvement. Either way, the remaining power cor-rections after subtracting the full NLL result shows a much steeper slope, which is as expected from its O(τ 2 cut ) scaling. This provides another example where considering only the overall size of the improvement can be potentially misleading. The gq + qg → Hq channel shows a similarly unusual behavior as for the leptonic definition.
In fig. 11 we show the rapidity dependence of the NLP corrections at fixed τ cut = 10 −3 for both leptonic and hadronic T normalized to the LO rapidity spectrum. The exponential enhancement for the hadronic definition at large |Y | is again apparent in the LL results. The NLL coefficients again exhibit an enhancement already for the leptonic definition at large Y . This is again due to the different x dependence of the quark PDF and the PDF derivatives compared to the LO gg luminosity to which we normalize. The quark PDF contributions are also the main reason why the NLL term for the gq channel (a 0 in table 2) is much larger than the LL contribution. For the hadronic definition, the e ±Y factors from the observable definition turn out to partially compensate these PDF effects. This is best visible in the gq channel, where the PDF enhanced terms at negative (LL) or positive (NLL) Y get reduced by a e ±Y factor from the observable definition. The same effect is also present in the gg channel at NLL. This is the reason why the a 0 term for the hadronic definition in the gq channel turns out to be even slightly smaller than for the leptonic definition.

Conclusions
In this paper, we have computed the next-to-leading power corrections in the N -jettiness resolution variable for Drell-Yan and gluon-fusion Higgs production at NLO. This builds on our previous work by computing the non-logarithmically enhanced terms at this order. These results enable the performance of the N -jettiness subtraction method to be improved, and provide important information on the structure of subleading power corrections beyond the leading logarithms. Our calculation is based on a master formula applicable to SCET I observables, and highlights a large degree of universality of these power corrections.
We explained in detail the issue of the treatment of Born measurements at subleading power. We have shown that an apparent disagreement in the literature arises due to the fact that the representation used to obtain the power corrections in refs. [59,66] is only valid when integrated over all rapidities, and therefore cannot be directly compared with the results of refs. [58,60] and those in the present paper, which are differential in rapidity. We show that after integration over rapidity the LL results agree. Further details can be found in sec. 6.
We find that the rapidity dependence of the NLL terms is quite sizeable and is therefore important to know to be able to improve the subtractions. One reason for this effect is the different x-dependence of parton luminosities or derivatives of PDFs appearing in the power corrections as compared to the Born-level parton luminosity. Hence, one can expect this to be a generic feature of subleading power corrections.
We also compared our analytic NLL results for gluon fusion Higgs production and Drell-Yan to numerical predictions for these NLO power corrections obtained from a fit to data from MCFM. In all cases, excellent agreement was found. In addition we studied the extent to which the inclusion of the NLL power corrections improves the subtraction. At NLO, the inclusion of the NLL power corrections completely captures the O(τ ) terms. Numerically, summing over production channels, the inclusion of these results reduces the size of the power corrections by two orders of magnitude in Higgs production and three orders of magnitude in Drell-Yan production.
There are a number of directions for future work. It will be interesting to extend the calculation of the NLL power corrections to NNLO. Generically we expect up to an order of magnitude improvement could also be obtained by extending the known LL power corrections at this order to NLL. Beyond fixed order, the derivation of subleading power renormalization group evolution equations at NLL would allow for the all-orders prediction of the NLL terms. Finally, while we have focused here on color-singlet production, our results provide an important step toward the calculation of the NLP corrections at higher orders and for more complicated processes. and S is the soft function. The beam functions can be further matched onto normal PDFs, All of these functions have definitions as field theory matrix elements in the EFT. Their fixed order definitions give rise to UV divergences, which are as usual removed by a renormalization procedure, which in turn gives rise to RGEs that can be used to resum large logarithms of T . In the approach presented in this paper, the same divergences appear as 1/ IR divergences in the soft and collinear limits of QCD amplitudes.
At LO, we have At one loop, the convolution structure thus becomes trivial. Working with hard, beam, soft, and PDFs in the bare factorization theorem we have Note the extra factor of Q in the beam contributions, arising from t a,b and T having different mass dimensions. We have written eq. (A.4) in a form similar to our master formulas, such that we can easily read off the one-loop beam function kernels and soft function. The arguments in eq. (A.4) all refer to ultraviolet divergences and can be removed by SCET counterterms to obtain the renormalized factorization theorem. To obtain this it is important to include virtual graphs in the various sectors as well as zero-bin subtractions for the beam functions.

A.1 Leading-Power Expansion of Matrix Elements
The leading-power behavior of real emission matrix elements in the soft and collinear limits is universal, see e.g. [7], and has already been used in sec. 4.4. Here, we briefly review the relevant formulas, and give the relevant one-loop expressions. Given the Born process κ a (q a ) + κ b (q b ) → L(q a + q b ) , (A.5) where the incoming momenta are given by q µ a = Qe Y n µ 2 , q µ b = Qe −Yn µ 2 , (A. 6) and we write the one-emission process as κ a (q a ) + κ b (q b ) → L(q a + q b − k) + κ 1 (k) . (A.7) In the soft limit k µ q µ a , q µ b , the squared matrix element obeys the LP relation where we made explicit that a soft emission can not change the incoming flavors. C = C F , C A is the Casimir constant for ab = qq, gg.
In the LP n-collinear limit, the particle k arises from the splitting κ a → κ a + κ 1 . If this splitting is allowed, at LP we have (in the notation of sec. 4) q a = q a /z a and q b = q b , and the LP limit of the matrix element is given by Similarly, in then-collinear limit arising from κ b → κ b + κ 1 , at LP we get q a = q a , q b = q b /z b , The one-loop splitting functions in d = 4 − 2 dimensions are given by [7] P qq (z, ) = C F 1 + z 2 1 − z − (1 − z) , P gq (z, ) = C F 1 + (1 − z) 2 z − z , Note that we flipped the notation of qg and gq relative to [7], following the standard convention.

A.2 NLO Soft Function
The NLO LP soft function follows from combining eq. (4.12) with eq. (A.8) using the same steps as in sec. 4.3, dσ (0,1) s The L n (x) are the standard one-dimensional plus distributions, see e.g. [84] for details. Note that there, the precise definition of the MS scheme is important. We use If one were to use µ 2 = (4π) Γ(1− ) µ 2 MS , one would miss the π 2 /3 term. For the NLP results presented in the main text, both definitions yield identical results.
Taking eq. (A.12) and adding the virtual soft diagram, and then comparing to eq. (A.4), the one-loop bare soft function can be read off as (A.14) The splitting function P ij (z) may contain divergences as z → 1, which are regulated by the overall (1 − z) − . All divergences thus arise from the two expansions where we defined µ 2 ρ = µ 2 ρ e Y for ease of notation. As written eq. (A. 17) does not yet contain the corresponding collinear virtual and zero-bin contributions.
Example: qq Kernel From eq. (A.11), we obtain where the LO quark splitting function is given by Adding the corresponding virtual collinear and zero-bin contributions, eq. (A.17) yields Note that all divergences are proportional to δ(1 − z), such that they cancel after adding the soft,n collinear and the virtual hard contribution from H NLO ab (Q 2 , ), as the latter also has the universal structure (for Drell-Yan) dσ (0,1) virt The cancellation of the 1/ 2 and the remaining 3/(2 ) pieces is obvious from comparing to eqs. (A.22) and (A.14). The P qq (z)/ term cancels with the ultraviolet divergence from the bare quark PDF. The remaining ln(Q 2 /µ 2 )/ term cancels when combining the L 0 (T /µ)/ and L 0 (t/µ 2 )/ distribution terms. The remaining O( 0 ) piece in eq. (A.22) gives the renormalized beam function and agrees with the result in ref. [38].